version: September 26, 2022

About this document


All analyses preformed with R version 4.1.0.

Basic setup of R environment


Loading required packages

For the following analyses we will require the use of a number of different R packages. We can use the following code to quickly load in the packages and install any packages not previously installed in the R console.

if (!require("pacman")) install.packages("pacman")

pacman::p_load_gh("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis", "ropensci/rnaturalearthhires")

pacman::p_load("cowplot", "car", "ggrepel", "ggspatial", "paletteer", "patchwork", "rgdal", "rnaturalearth", "sf", "Hmisc", "MCMC.OTU", "pairwiseAdonis", "RColorBrewer", "Redmonder", "flextable", "lubridate", "officer", "adegenet", "dendextend", "gdata", "ggdendro", "hierfstat", "Imap", "kableExtra", "poppr", "reshape2", "StAMPP", "vcfR", "vegan", "boa", "magick", "rgeos", "sdmpredictors", "ggcorrplot", "tidyverse", "TeachingDemos", "LaplacesDemon", "adespatial", "ggnewscale")

options("scipen" = 10)

load("fknmsSint.RData")
# save.image("fknmsSint.RData")


Sampling info


Map of study sites


fknmsSites = read.csv("../data/stephanocoeniaMetaData.csv", header = TRUE)
fknmsSites$depthZone = factor(fknmsSites$depthZone)
fknmsSites$depthZone = factor(fknmsSites$depthZone, levels = levels(fknmsSites$depthZone)[c(2,1)])

fknmsSites$site = factor(fknmsSites$site)
fknmsSites$site = factor(fknmsSites$site, levels = levels(fknmsSites$site)[c(4, 1, 3, 2)])
fknmsSites$date = mdy(fknmsSites$date) %>% format("%d %b %Y")

fknmsPops = fknmsSites %>% group_by(site) %>% summarise(latDD = mean(latDD), longDD = mean(longDD), n = n()) %>% droplevels()

fknmsSampleSites = fknmsSites %>% group_by(site, siteID, depthZone) %>% summarise(latDD = min(latDD), longDD = min(longDD))
## `summarise()` has grouped output by 'site', 'siteID'. You can override using the
## `.groups` argument.
fknmsBounds = read.csv("../data/fknmsSPA.csv", header = TRUE)

states = st_as_sf(ne_states(country = c("United States of America")), scale = "count",  crs = 4326) %>% filter(name_en %in% c("Florida", "Georgia", "Alabama"))
countries = st_as_sf(ne_countries(country = c("Cuba", "Mexico", "The Bahamas", "Bermuda"), scale = "Large"), crs = 4326)
bahamas = read_sf("../data/shp/bahamasShoreline.shp") %>% st_transform(crs = 4326)
cuba = read_sf("../data/shp/cubaShoreline.shp") %>% st_transform(crs = 4326)
florida = read_sf("../data/shp/floridaShoreline.shp") %>% st_transform(crs = 4326)
bathy = read_sf("../data/shp/flBathy.shp") %>% st_transform(crs = 4326) %>% subset(subset = DATASET %in% c("fl_shelf", "fl_coast"))
tortugasBathy = read_sf("../data/shp/tortugasBathy.shp") %>% st_transform(crs = 4326)


Next we build a hi-res polygon of FL with the study site marked and a zoomed in map of the colony locations. We use ggspatial to add a north arrow and scale bar to the main map.

flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]
boundPal = c("gray30", paletteer_d("vapoRwave::vapoRwave")[10])
pink = paletteer_d("vapoRwave::vapoRwave")[5]
purple = paletteer_d("vapoRwave::vapoRwave")[10]
  
floridaMap = ggplot() +
  geom_polypath(data = fknmsBounds[fknmsBounds$type == "Sanctuary",], aes(x = long, y = lat, group = location), alpha = 0.1, fill = "black", color = NA) +
  geom_polypath(data = fknmsBounds[fknmsBounds$location == "FKNMS2",], aes(x = long, y = lat), fill = "aliceblue", color = NA) +
  geom_polypath(data = fknmsBounds, aes(x = long, y = lat, color = type, group = location), fill = NA) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = boundPal, name = "Boundaries", labels = c("FKNMS", "SPA")) +
  geom_point(data = fknmsSites, aes(x = longDD, y = latDD, shape = depthZone), size = 0) +
  scale_shape_manual(values = c(24, 25), name = "Depth") +
  geom_sf(data = florida, fill = "white", size = 0.15) +
  geom_sf(data = cuba, fill = "white", size = 0.15) +
  geom_sf(data = bahamas, fill = "white", size = 0.15) +
  geom_segment(aes(x = -80.1, y = 25.3, xend = -78.825, yend = 24.45), size = 0.25) +
  geom_segment(aes(x = -80.4, y = 25, xend = -80.27, yend = 23), size = 0.25) +
  geom_segment(aes(x = -81.75, y = 24.7, xend = -82.22, yend = 24.29), size = 0.25) +
  geom_segment(aes(x = -81.45, y = 24.7, xend = -80.77, yend = 24.29), size = 0.25) +
  geom_segment(aes(x = -83.25, y = 24.75, xend = -84.183, yend = 24.29), size = 0.25) +
  geom_segment(aes(x = -82.95, y = 24.75, xend = -82.73, yend = 24.29), size = 0.25) +
  geom_rect(aes(xmin = -80.4, xmax = -80.1, ymin = 25, ymax = 25.3), fill = pink, color = "black", size = 0.35, alpha = 0.5) +
  geom_rect(aes(xmin = -81.75, xmax = -81.45, ymin = 24.4, ymax = 24.7), fill = pink, color = "black", size = 0.25, alpha = 0.5) +
  geom_rect(aes(xmin = -83.25, xmax = -82.95, ymin = 24.45, ymax = 24.75), fill = pink, color = "black", size = 0.25, alpha = 0.5) +
  geom_point(data = fknmsPops, aes(x = longDD, y = latDD, fill = site), shape = 21, size = 2) +
  coord_sf(xlim = c(-84, -79), ylim = c(23, 27)) +
  scale_x_continuous(breaks = c(seq(-84, -79, by = 1))) +
  scale_y_continuous(breaks = c(seq(23, 27, by = 1))) +
  annotation_scale(location = "br", pad_x = unit(1.35, "cm"), text_pad = unit(-3.75, "cm")) +
  guides(fill = guide_legend(override.aes = list(shape = 22, color = NA, size = 4), ncol = 2, order = 1), shape = guide_legend(override.aes = list(size = 3), order = 2), color = guide_legend(override.aes = list(fill = "black", alpha = 0.1), order = 3)) +
  theme_bw() +
  theme(panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        axis.title = element_blank(),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.box = "horizontal")

# floridaMap

largeMap = inset = ggplot() +
  geom_sf(data = states, fill = "white", size = 0.3) +
  geom_sf(data = countries, fill = "white", size = 0.3) +
  geom_rect(aes(xmin = -84, xmax = -79, ymin = 23, ymax = 27), color = pink, fill = NA, size = 0.75) +
  geom_rect(aes(xmin = -78.8, xmax = -77, ymin = 22.2, ymax = 22.6), fill = "aliceblue", color = NA) +
  annotation_scale(location = "bl", pad_x = unit(1.85, "cm")) +
  annotation_north_arrow(location = "tr", style = north_arrow_minimal(), pad_x = unit(-0.3, "cm")) +
  coord_sf(xlim = c(-87, -76), ylim = c(22, 31)) +
  theme_bw() +
  theme(legend.title = element_text(size = 9, face = "bold"),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        legend.text = element_text(size = 9),
        legend.position = "none",
        plot.background = element_blank())

# largeMap

inset = ggplot() +
  geom_polypath(data = fknmsBounds[fknmsBounds$type == "Sanctuary",], aes(x = long, y = lat, group = location), alpha = 0.1, fill = "black", color = NA) +
   geom_polypath(data = fknmsBounds[fknmsBounds$location == "FKNMS2",], aes(x = long, y = lat), fill = "aliceblue", color = NA) +
  geom_segment(aes(x = -82.9645, xend = -82.4, y = 24.6, yend = 24.6), color = "gray92", size = .55) +
  geom_sf(data = bathy, color = "gray75", size = 0.25) +
  geom_polypath(data = fknmsBounds, aes(x = long, y = lat, color = type, group = location), fill = NA) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = boundPal, name = "Boundaries", labels = c("FKNMS", "SPA")) +
  geom_point(data = fknmsSampleSites, aes(x = longDD, y = latDD, fill = site, shape = depthZone), size = 1.5) +
  geom_sf(data = florida, fill = "white", size = 0.15) +
  scale_shape_manual(values = c(24, 25), name = "Depth") +
  theme_bw() +
  theme(legend.title = element_text(size = 9, face = "bold"),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        legend.text = element_text(size = 9),
        legend.position = "none",
        plot.background = element_blank())

# inset

upperKeys = inset +
  # annotation_scale(location = "br") +
  annotation_scale(location = "bl", pad_x = unit(1.9, "cm")) +
  coord_sf(xlim = c(-80.4, -80.1), ylim = c(25.0, 25.3)) +
  scale_x_continuous(breaks = c(seq(-80.4, -80.0, by = .1))) +
  scale_y_continuous(breaks = c(seq(25.0, 25.3, by = .1)))

lowerKeys = inset +
  # annotation_scale(location = "br") +
  annotation_scale(location = "bl", pad_x = unit(1.9, "cm")) +
  coord_sf(xlim = c(-81.75, -81.45), ylim = c(24.4, 24.7)) +
  scale_x_continuous(breaks = c(seq(-81.7, -81.3, by = .1))) +
  scale_y_continuous(breaks = c(seq(24.4, 24.7, by = .1)))

dryTortugas = ggplot() +
  geom_polypath(data = fknmsBounds[fknmsBounds$type == "Sanctuary",], aes(x = long, y = lat, group = location), alpha = 0.1, fill = "black", color = NA) +
   geom_polypath(data = fknmsBounds[fknmsBounds$location == "FKNMS2",], aes(x = long, y = lat), fill = "aliceblue", color = NA) +
  geom_segment(aes(x = -82.9645, xend = -82.4, y = 24.6, yend = 24.6), color = "gray92", size = .55) +
  geom_sf(data = tortugasBathy, color = "gray75", size = 0.25) +
  geom_polypath(data = fknmsBounds, aes(x = long, y = lat, color = type, group = location), fill = NA) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = boundPal, name = "Boundaries", labels = c("FKNMS", "SPA")) +
  geom_point(data = fknmsSites, aes(x = longDD, y = latDD, fill = site, shape = depthZone), size = 1.5) +
  geom_sf(data = florida, fill = "white", size = 0.25) +
  scale_shape_manual(values = c(24, 25), name = "Depth") +
  # annotation_scale(location = "br", pad_x = unit(1.9, "cm")) +
  annotation_scale(location = "bl", pad_x = unit(1.9, "cm")) +
  # coord_sf(xlim = c(-83.25, -82.95), ylim = c(24.45, 24.75)) +
  coord_sf(xlim = c(-83.25, -82.95), ylim = c(24.45, 24.75)) +
  scale_x_continuous(breaks = c(seq(-83.2, -82.9, by = .1))) +
  scale_y_continuous(breaks = c(seq(24.4, 24.7, by = .1))) +
  theme_bw() +
  theme(legend.title = element_text(size = 9, face = "bold"),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        legend.text = element_text(size = 9),
        legend.position = "none",
        plot.background = element_blank())

stephanoPic = image_read("../data/stephano.png") %>%
  image_border("black", "23x23")

fknmsMap = ggdraw() +
  draw_plot(floridaMap) +
  draw_plot(upperKeys, x = 0.695, y = 0.22, width = 0.29, height = 0.29) +
  draw_plot(lowerKeys, x = 0.38, y = 0.19, width = 0.29, height = 0.29) +
  draw_plot(dryTortugas, x = 0.062, y = 0.19, width = 0.29, height = 0.29) +
  # draw_plot(largeMap, x = 0.685, y = 0.712, width = 0.29, height = 0.29) +
  # draw_image(stephanoPic, x = 0.092, y = 0.732, width = 0.24, height = 0.254)
  draw_plot(largeMap, x = 0.075, y = 0.712, width = 0.29, height = 0.29) +
  draw_image(stephanoPic, x = 0.722, y = 0.732, width = 0.24, height = 0.254)

# fknmsMap

ggsave("../figures/figure1.png", plot = fknmsMap, width = 16, height = 16, units = "cm", dpi = 300)
ggsave("../figures/figure1.svg", plot = fknmsMap, width = 16, height = 16, units = "cm", dpi = 300)


S. intersepta population genetics from SNPs


Analyzing 2bRAD generated SNPs (24,670 loci) for population structure//genetic connectivity across sites and depth zones in FKNMS

How many reads?

rawSintReads = read.delim("../data/snps/sintRawReadCounts", header = FALSE)
colnames(rawSintReads) = c("sample", "reads")

head(rawSintReads)
##    sample    reads
## 1 FKSi1-1 42167284
## 2 FKSi1-2 54651139
## 3 FKSi1-3 41635251
## 4 FKSi1-4 37754282
## 5 FKSi1-5 39973126
## 6 FKSi1-6 45580831
#total reads
sum(rawSintReads$reads)
## [1] 796139328
#average reads/sample
(sum(rawSintReads$reads)/226)
## [1] 3522740

Identifiying clonal multi-locus genotypes

Dendrogram with clones

Identification of any natural clones using technical replicates as a baseline for clonality between samples.

cloneBams = read.csv("../data/stephanocoeniaMetaData.csv") # list of bam files

cloneMa = as.matrix(read.table("../data/snps/clones/sintClones.ibsMat")) # reads in IBS matrix produced by ANGSD 

dimnames(cloneMa) = list(cloneBams[,1],cloneBams[,1])
# clonesHc = hclust(as.dist(cloneMa),"ave")

clonePops = cloneBams$site
cloneDepth = cloneBams$depthZone

cloneDend = cloneMa %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
cloneDData = cloneDend %>% dendro_data()

# Making the branches hang shorter so we can easily see clonal groups
cloneDData$segments$yend2 = cloneDData$segments$yend
for(i in 1:nrow(cloneDData$segments)) {
  if (cloneDData$segments$yend2[i] == 0) {
    cloneDData$segments$yend2[i] = (cloneDData$segments$y[i] - 0.01)}}

cloneDendPoints = cloneDData$labels
cloneDendPoints$pop = clonePops[order.dendrogram(cloneDend)]
cloneDendPoints$depth=cloneDepth[order.dendrogram(cloneDend)]
rownames(cloneDendPoints) = cloneDendPoints$label

# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(cloneDData$segments)) {
  if (cloneDData$segments$yend[i] == 0) {
    point[i] = cloneDData$segments$y[i] - 0.01
  } else {
    point[i] = NA}}

cloneDendPoints$y = point[!is.na(point)]

techReps = c("SFK066.1", "SFK066.2", "SFK066.3", "SFK162.1", "SFK162.2", "SFK162.3", "SFK205.1", "SFK205.2", "SFK205.3")

cloneDendPoints$depth = factor(cloneDendPoints$depth)
cloneDendPoints$depth = factor(cloneDendPoints$depth, levels(cloneDendPoints$depth)[c(2,1)])

cloneDendPoints$pop = factor(cloneDendPoints$pop)
cloneDendPoints$pop = factor(cloneDendPoints$pop,levels(cloneDendPoints$pop)[c(4, 1, 3, 2)])

flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]

cloneDendA = ggplot() +
  geom_rect(aes(xmin = 44.25, xmax = 47.75, ymin = 0.03, ymax = 0.085), fill = pink, alpha = 0.4) +
  geom_rect(aes(xmin = 154.25, xmax = 157.75, ymin = 0.065, ymax = 0.12), fill = pink, alpha = 0.4) +
  geom_rect(aes(xmin = 172.25, xmax = 175.75, ymin = 0.065, ymax = 0.12), fill = pink, alpha = 0.4) +
  geom_segment(data = segment(cloneDData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
  geom_point(data = cloneDendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
  scale_fill_manual(values = flPal, name= "Site:") +
  scale_shape_manual(values = c(24, 25), name = "Depth Zone:") +
  geom_hline(yintercept = 0.118, color = pink, lty = 5, size = 1) + # creating a dashed line to indicate a clonal distance threshold
  geom_text(data = subset(cloneDendPoints, subset = label %in% techReps), aes(x = x, y = (y - .02), label = label), angle = 90) + # spacing technical replicates further from leaf
  geom_text(data = subset(cloneDendPoints, subset = !label %in% techReps), aes(x = x, y = (y - .015), label = label), angle = 90) +
  labs(y = "Genetic distance (1 - IBS)") +
  guides(fill = guide_legend(override.aes = list(shape = 22, size = 10), ncol = 2), shape = guide_legend(override.aes = list(size = 8), ncol = 1, order = 1)) +
  coord_cartesian(xlim = c(5, 218), ylim = c(0.03, 0.25)) +
  theme_classic()

cloneDend = cloneDendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 24, color = "black", angle = 90),
  axis.text.y = element_text(size = 20, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 24),
  legend.text = element_text(size = 20),
  legend.position = "bottom")

# cloneDend

# ggsave("../figures/rmd/cloneDend.png", plot = cloneDend, height = 8, width = 35, units = "in", dpi = 300)


We removed the technical replicates/clones and re-ran ANGSD for all the following pop-gen analyses.

Dendrogram without clones

bams = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66,68,164,166,209,211),] # list of bams files and their populations (tech reps removed)

ma = as.matrix(read.table("../data/snps/sintNoClones.ibsMat")) # reads in IBS matrix produced by ANGSD 

dimnames(ma) = list(bams[,1],bams[,1])

pops = bams$site
depth = bams$depthZone

dend = ma %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
dData = dend %>% dendro_data()

# Making the branches hang shorter so we can easily see clonal groups
dData$segments$yend2 = dData$segments$yend
for(i in 1:nrow(dData$segments)) {
  if (dData$segments$yend2[i] == 0) {
    dData$segments$yend2[i] = (dData$segments$y[i] - 0.01)}}

dendPoints = dData$labels
dendPoints$pop = pops[order.dendrogram(dend)]
dendPoints$depth = depth[order.dendrogram(dend)]
rownames(dendPoints) = dendPoints$label

# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(dData$segments)) {
  if (dData$segments$yend[i] == 0) {
    point[i] = dData$segments$y[i] - 0.01
  } else {
    point[i] = NA}}

dendPoints$y = point[!is.na(point)]

dendPoints$depth = factor(dendPoints$depth)
dendPoints$depth = factor(dendPoints$depth, levels(dendPoints$depth)[c(2,1)])

dendPoints$pop = factor(dendPoints$pop)
dendPoints$pop = factor(dendPoints$pop, levels(dendPoints$pop)[c(4, 1, 3, 2)])

flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]

dendA = ggplot() +
  geom_segment(data = segment(dData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
  geom_point(data = dendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
  geom_text(data = dendPoints, aes(x = x, y = (y - .01), label = label), angle = 90) + 
  scale_fill_manual(values = flPal, name= "Site:")+
  scale_shape_manual(values = c(24, 25), name = "Depth Zone:")+
 # spacing technical replicates further from leaf
  labs(y = "Genetic distance (1 - IBS)") +
  guides(fill = guide_legend(override.aes = list(shape = 22, size = 10), ncol = 2), shape = guide_legend(override.aes = list(size = 8), ncol = 1, order = 1)) +
  coord_cartesian(xlim = c(5, 218)) +
  theme_classic()

dend = dendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 24, color = "black", angle = 90),
  axis.text.y = element_text(size = 20, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 24),
  legend.text = element_text(size = 20),
  legend.position = "bottom")

# dend

# ggsave("../figures/rmd/dend.png", plot = dend, height = 6, width = 37, units = "in", dpi = 300)


Dendrogram plots

dendPlots = (cloneDend / dend) +
  plot_annotation(tag_levels = 'A') +
  plot_layout(guides = "collect") & 
  theme(plot.tag = element_text(size = 32),
        legend.position = "bottom")

ggsave("../figures/figureS1.png", plot = dendPlots, height = 13, width = 35, units = "in", dpi = 300)
ggsave("../figures/figureS1.svg", plot = dendPlots, height = 13, width = 35, units = "in", dpi = 300)


Population statistics

Heterozygosity and Inbreeding

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "Site" = site, "Depth" = depthZone) # Reads in population data
popData$a = c(0:219)

popData$Site = factor(popData$Site)
popData$Site = factor(popData$Site, levels = levels(popData$Site)[c(2,3,1,4)]) 
popData$Depth = factor(popData$Depth)
popData$Depth = factor(popData$Depth, levels = levels(popData$Depth)[c(2,1)]) 

sampleData = fknmsSites[-c(66,68,164,166,209,211),] %>% group_by(site, depthZone)%>% summarise(depthZone = (first(depthZone)), depthRange = paste(min(depthM), "--", max(depthM), sep = ""), meanDepth = round(mean(depthM),1), n = n())%>% droplevels() %>% as.data.frame()
## `summarise()` has grouped output by 'site'. You can override using the `.groups`
## argument.
# Average depth of populations
fkPopDepths = fknmsSites[-c(66,68,164,166,209,211),] %>%  group_by(site, depthZone) %>% summarise(avgDepthM = mean(depthM), n = n())
## `summarise()` has grouped output by 'site'. You can override using the `.groups`
## argument.
fkPopDepths
## # A tibble: 8 × 4
## # Groups:   site [4]
##   site          depthZone  avgDepthM     n
##   <fct>         <fct>          <dbl> <int>
## 1 Upper Keys    Shallow         23.6    30
## 2 Upper Keys    Mesophotic      43.8    30
## 3 Lower Keys    Shallow         18.0    30
## 4 Lower Keys    Mesophotic      32.8    30
## 5 Tortugas Bank Shallow         21.1    30
## 6 Tortugas Bank Mesophotic      32.0    25
## 7 Riley's Hump  Shallow         26.4    30
## 8 Riley's Hump  Mesophotic      33.2    15
sampleTab = sampleData
colnames(sampleTab) = c("Site", "Depth zone", "Sampling \ndepth (m)", "Average \ndepth (m)", "n")

sampleTab$Site
## [1] Upper Keys    Upper Keys    Lower Keys    Lower Keys    Tortugas Bank Tortugas Bank
## [7] Riley's Hump  Riley's Hump 
## Levels: Upper Keys Lower Keys Tortugas Bank Riley's Hump
finalTabSite = c("Upper Keys", "", "Lower Keys","", "Tortugas Bank", "", "Riley's Hump", "")

sampleTab$Site = finalTabSite

hetAll = read.table("../data/snps/hetAllSites") 
colnames(hetAll) = c("sample", "All")
hetAll$sample = str_pad(hetAll$sample, 3, pad = "0")
hetAll$sample = paste("SFK",hetAll$sample, sep ="")

hetSnps = read.table("../data/snps/hetSnps")
colnames(hetSnps) = c("sample", "SNPs")
hetSnps$sample = str_pad(hetSnps$sample, 3, pad = "0")
hetSnps$sample = paste("SFK",hetSnps$sample, sep ="")

sintBreed = read.delim("../data/snps/newres")
  
sintBreed2 = sintBreed %>% group_by(a) %>% dplyr::select("inbreed" = Fa)
## Adding missing grouping variables: `a`
sintBreed3 = sintBreed %>% group_by(b) %>% dplyr::select("inbreed" = Fb)
## Adding missing grouping variables: `b`
sintBreed = bind_rows(sintBreed2, sintBreed3) %>% group_by(a) %>% summarise("inbreed" = mean(inbreed)) 
sintRelate = read.delim("../data/snps/newres")
sintRelate2 = sintRelate %>% group_by(a, b) %>% dplyr::select("Rab" = rab, "theta" = theta)
## Adding missing grouping variables: `a`, `b`
sintRelate2 = sintRelate2 %>% left_join(popData, by = "a") %>% left_join(popData, by = c("b" = "a"), suffix = c(".a", ".b")) 

sintRelate2$popDepth.a = paste(sintRelate2$Site.a, sintRelate2$Depth.a, sep = " ")
sintRelate2$popDepth.b = paste(sintRelate2$Site.b, sintRelate2$Depth.b, sep = " ")

sintRelate = sintRelate2 %>% filter(popDepth.a == popDepth.b) %>% rename(Depth = Depth.a, Site = Site.a)

sintRelateMean =  sintRelate %>% group_by(Site, Depth) %>% dplyr::summarize(N = n(), meanRab = mean(Rab), seRab = sd(Rab)/sqrt(N), meanTheta = mean(theta), seTheta = sd(theta)/sqrt(N)) %>% dplyr::select(-N)
## `summarise()` has grouped output by 'Site'. You can override using the `.groups`
## argument.
het = left_join(popData, hetAll, by = "sample") %>% left_join(hetSnps, by = "sample") %>% mutate("inbreed" = sintBreed$inbreed) 

hetStats = het %>% group_by(Site, Depth) %>% summarise(N = n(), meanAll = mean(All), sdAll = sd(All), seAll = sd(All)/sqrt(N), meanSnps = mean(SNPs), sdSnps = sd(SNPs), seSnps = sd(SNPs)/sqrt(N), meanInbreed = mean(inbreed), sdInbreed = sd(inbreed), seInbreed = sd(inbreed)/sqrt(N)) %>% left_join(sintRelateMean)
## `summarise()` has grouped output by 'Site'. You can override using the `.groups`
## argument.
## Joining, by = c("Site", "Depth")
min(hetStats$meanAll, na.rm = TRUE)
## [1] 0.00406
max(hetStats$meanAll, na.rm = TRUE)
## [1] 0.004352
min(hetStats$meanSnps, na.rm = TRUE)
## [1] 0.1999967
max(hetStats$meanSnps, na.rm = TRUE)
## [1] 0.23806
hetTab = hetStats %>% arrange(desc(Site))

hetTab$n = hetTab$N
hetTab$Ha = paste(round(hetTab$meanAll, 4), "±", round(hetTab$seAll, 5), sep = " ")
hetTab$Hv = paste(round(hetTab$meanSnps, 3), "±", round(hetTab$seSnps, 4), sep = " ")
hetTab$F =  paste(round(hetTab$meanInbreed, 2), "±", round(hetTab$seInbreed, 3), sep = " ")
hetTab$Rab =  paste(round(hetTab$meanRab, 2), "±", round(hetTab$seRab, 3), sep = " ")
hetTab$Theta =  paste(round(hetTab$meanTheta, 2), "±", round(hetTab$seTheta, 4), sep = " ")

hetTab$`Sampling \ndepth (m)` = sampleTab$`Sampling \ndepth (m)`
hetTab$`Average \ndepth (m)` = sampleTab$`Average \ndepth (m)`

hetTab = hetTab %>% dplyr::select(Site, Depth, `Sampling \ndepth (m)`, `Average \ndepth (m)`, n, Ha, Hv, F, Rab, Theta)
colnames(hetTab)[2] = "Depth \nzone"

finalTabSite = c("Upper Keys", "", "Lower Keys", "", "Tortugas Bank", "", "Riley's Hump", "")

hetTab$Site = finalTabSite

hetTabPub = hetTab %>% dplyr::select(-Theta) %>% 
  flextable() %>%
  flextable::compose(part = "header", j = "n", value = as_paragraph(as_i("n"))) %>%
   flextable::compose(part = "header", j = "Ha", value = as_paragraph(as_i("H"), as_sub("A"))) %>%
  flextable::compose(part = "header", j = "Hv", value = as_paragraph(as_i("H"),as_sub("V"))) %>%
  flextable::compose(part = "header", j = "F", value = as_paragraph(as_i("F"))) %>%
  # flextable::compose(part = "header", j = "Rab", value = as_paragraph(as_i("R"), as_i(as_sub("AB")))) %>%
  flextable::compose(part = "header", j = "Rab", value = as_paragraph(as_i("R"), as_i(as_sub("AB")))) %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::fontsize(size = 10, part = "all") %>%
  flextable::bold(part = "header") %>%
  flextable::align(align = "left", part = "all") %>%
  flextable::autofit()

table2 = read_docx()
table2 = body_add_flextable(table2, value = hetTabPub)
print(table2, target = "../tables/table2.docx")

hetTabPub

ANOVA on population statistics

hetAllLm = lm(data = het, All~Site*Depth)
hetSnpLm = lm(data = het, SNPs~Site*Depth)
inbreedLm = lm(data = het, inbreed~Site*Depth)
relateLm = lm(data = sintRelate, Rab~Site*Depth)
kinLm = lm(data = sintRelate, theta~Site*Depth)

hetAllANOVA = summary(aov(hetAllLm))
hetAllANOVA
##              Df      Sum Sq      Mean Sq F value  Pr(>F)   
## Site          3 0.000000785 0.0000002618   1.833 0.14210   
## Depth         1 0.000001133 0.0000011330   7.936 0.00531 **
## Site:Depth    3 0.000000103 0.0000000344   0.241 0.86760   
## Residuals   212 0.000030268 0.0000001428                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hetSnpANOVA = summary(aov(hetSnpLm))
hetSnpANOVA
##              Df  Sum Sq Mean Sq F value         Pr(>F)    
## Site          3 0.00115 0.00038   0.370          0.775    
## Depth         1 0.04412 0.04412  42.514 0.000000000505 ***
## Site:Depth    3 0.00355 0.00118   1.139          0.334    
## Residuals   212 0.22002 0.00104                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
inbreedANOVA = summary(aov(inbreedLm))
inbreedANOVA
##              Df Sum Sq Mean Sq F value            Pr(>F)    
## Site          3 0.0115  0.0038   0.460             0.711    
## Depth         1 0.5135  0.5135  61.719 0.000000000000196 ***
## Site:Depth    3 0.0395  0.0132   1.581             0.195    
## Residuals   212 1.7637  0.0083                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
relateANOVA = summary(aov(relateLm))
relateANOVA
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Site           3   0.06  0.0196   1.764    0.152    
## Depth          1   0.75  0.7537  67.736 2.75e-16 ***
## Site:Depth     3   0.05  0.0171   1.538    0.202    
## Residuals   3007  33.46  0.0111                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Individual relatedness heatmap

sintFam = sintRelate2 
sintFam$popDepth.a = factor(sintFam$popDepth.a)
sintFam$popDepth.b = factor(sintFam$popDepth.b)

sintFam$popDepth.a = factor(sintFam$popDepth.a, levels = levels(sintFam$popDepth.a)[c(8, 7, 2, 1, 6, 5, 4, 3)])
sintFam$popDepth.b = factor(sintFam$popDepth.b, levels = levels(sintFam$popDepth.b)[c(8, 7, 2, 1, 6, 5, 4, 3)])

sintMatDat = sintFam %>% as.data.frame () %>% dplyr::select(sample.a, sample.b, Rab)
nameVals <- sort(unique(unlist(sintMatDat[1:2])))
# construct 0 matrix of correct dimensions with row and column names
sintMat <- matrix(0, length(nameVals), length(nameVals), dimnames = list(nameVals, nameVals))

# fill in the matrix with matrix indexing on row and column names
sintMat[as.matrix(sintMatDat[c("sample.a", "sample.b")])] <- sintMatDat[["Rab"]]
lowerTriangle(sintMat, byrow = TRUE) <- upperTriangle(sintMat)

sintFam = sintFam %>% arrange(popDepth.a)

sintFam$sample.a = factor(sintFam$sample.a)
sintFam$sample.a = factor(sintFam$sample.a, levels = unique(sintFam$sample.a))
sintFam$sample.a = factor(sintFam$sample.a, levels = c("SFK220", levels(sintFam$sample.a)))
sintFam$sample.b = factor(sintFam$sample.b)
sintFam$sample.b = factor(sintFam$sample.b, levels = levels(sintFam$sample.a))

sampleOrder = levels(sintFam$sample.a)

sintMat <- sintMat[,sampleOrder] %>%
  .[sampleOrder,]

diag(sintMat)<-NA

sintMat <- as.data.frame(sintMat)

sintMat$sample.a = factor(colnames(sintMat))

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample.a" = tubeID, "site" = site, "depth" = depthZone) %>% mutate(popDepth = paste(site, depth, sep = "\\n"))

sintMat = sintMat %>% left_join(popData)
## Joining, by = "sample.a"
sintMat$sample.a = factor(sintMat$sample.a, levels = unique(sampleOrder))

sintFam = melt(sintMat, id.vars = c("sample.a", "site", "depth", "popDepth"), value.name = "relate", variable.name = "sample.b", na.rm = FALSE)

sintFam$site = factor(sintFam$site)
sintFam$site = factor(sintFam$site, levels = levels(sintFam$site)[c(4, 1, 3, 2)])
sintFam$depth = factor(sintFam$depth)
sintFam$depth = factor(sintFam$depth, levels = levels(sintFam$depth)[c(2, 1)])

relateHeatmapA = ggplot(data = sintFam, aes(sample.a, sample.b, fill = relate)) +
  geom_tile(color = "white") +
  geom_segment(data = sintFam, aes(x = sample.a, xend = sample.a, y = 0, yend = -4, color = site), size = 1) +
  geom_segment(data = sintFam, aes(x = 0, xend = -4, y = sample.a, yend = sample.a, color = site), size = 1) +
  scale_color_manual(values = flPal, name = "Site", guide = guide_legend(override.aes = list(size = 7), order = 2)) +
  geom_vline(xintercept = c(30.5, 60.5, 90.5, 120.5, 150.5, 175.5, 205.5), size = 0.1, color = "black") +
  geom_hline(yintercept = c(15.5, 45.5, 70.5, 100.5, 130.5, 160.5, 190.5), size = 0.1, color = "black") +
  new_scale("color") +
  geom_segment(data = sintFam, aes(x = sample.a, xend = sample.a, y = -2, yend = -4, color = depth), size = 1) +
  geom_segment(data = sintFam, aes(x = -2, xend = -4, y = sample.a, yend = sample.a, color =      depth), size = 1) +
  scale_color_manual(values = c("gray75", "gray35"), name = "Depth Zone", guide = guide_legend(override.aes = list(size = 7), order = 3)) +
  scale_fill_gradientn(colors = paletteer_d("RColorBrewer::RdPu"), space = "Lab", name = "Relatedness", na.value = "black", limits = c(0,0.6), breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6), guide = guide_colorbar(order = 1)) +
  scale_y_discrete(limits = rev(levels(sintFam$sample.b))) +
  coord_cartesian(xlim = c(-3.5,220)) +
  theme_minimal()

relateHeatmap = relateHeatmapA + theme(
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  panel.grid.major = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  legend.title = element_text(size = 14),
  legend.text = element_text(size = 14),
  plot.title = element_text(size = 16)
)

# relateHeatmap

ggsave("../figures/relateHeatmap.png", plot = relateHeatmap, width = 27, height = 22, units = "cm", dpi = 300)
ggsave("../figures/relateHeatmap.svg", plot = relateHeatmap, width = 27, height = 22, units = "cm", dpi = 300)


Relatedness dendrogram

sintFam2 = sintRelate2 %>% filter(sample.a != "NA") %>% dplyr::select("sample.a" = a, "sample.b" = b, "relate" = Rab) %>% as.data.frame () %>% mutate(sample.a = paste("SFK", str_pad(sample.a +1, 3, pad = "0" ), sep = ""), sample.b = paste("SFK", str_pad(sample.b +1, 3, pad = "0" ), sep = ""))

# sintFam2$relate = sintFam2$relate %>% replace_na(1)
# sintFam2$relate = 1-sintFam2$relate

nameVals <- sort(unique(unlist(sintFam2[1:2])))
# construct 0 matrix of correct dimensions with row and column names
relateMa <- matrix(0, length(nameVals), length(nameVals), dimnames = list(nameVals, nameVals))

relateMa[as.matrix(sintFam2[c("sample.a", "sample.b")])] <- sintFam2[["relate"]]
lowerTriangle(relateMa, byrow = TRUE, diag = FALSE) <- upperTriangle(relateMa, diag = FALSE)
diag(relateMa) = 1

relateDend = hclust(as.dist(1-relateMa), "single") %>% as.dendrogram()
relateDData = relateDend %>% dendro_data()

# Making the branches hang shorter so we can easily see clonal groups
relateDData$segments$yend2 = relateDData$segments$yend
for(i in 1:nrow(relateDData$segments)) {
  if (relateDData$segments$yend2[i] == 0) {
    relateDData$segments$yend2[i] = (relateDData$segments$y[i] - 0.05)}}

relBams = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66,68,164,166,209,211),] %>% dplyr::select("label" = Sample, "pop" = site, "depth" = depthZone)
relBams$label = gsub(pattern = "\\.2", "", relBams$label)

relDendPoints = relateDData$labels %>% left_join(relBams)
## Joining, by = "label"
rownames(relDendPoints) = relDendPoints$label

relDendPoints$depth = factor(relDendPoints$depth)
relDendPoints$depth = factor(relDendPoints$depth, levels(relDendPoints$depth)[c(2,1)])

relDendPoints$pop = factor(relDendPoints$pop)
relDendPoints$pop = factor(relDendPoints$pop, levels(relDendPoints$pop)[c(4, 1, 3, 2)])

# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(relateDData$segments)) {
  if (relateDData$segments$yend[i] == 0) {
    point[i] = relateDData$segments$y[i] - 0.05
  } else {
    point[i] = NA}}

relDendPoints$y = point[!is.na(point)]

#Assign dominant lineage (75% or greater)
domClust = fkSintAdmix %>% dplyr::select(-"cluster2.1", -"cluster2.2", -"ord", -"pop", -"depth", -"popdepth") %>% left_join(het) %>% mutate(cluster = ifelse(cluster4.1 < 0.75 & cluster4.2 < 0.75  & cluster4.3 < 0.75 & cluster4.4 < 0.75, "NA", ifelse(cluster4.1 >=0.75, 1, ifelse(cluster4.2 >= 0.75, 2, ifelse(cluster4.3 >= 0.75, 3,ifelse(cluster4.4 >= 0.75, 4, 0))))))
## Joining, by = "sample"
domClust$cluster = as.factor(domClust$cluster)
levels(domClust$cluster) = c("Red", "Brown", "Tan", "Orange", "Admixed")

relDendPoints2 = domClust %>% dplyr::select("label" = sample, cluster) %>% left_join(relDendPoints)
## Joining, by = "label"
flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]
kColPal = c("#A14747", "#A3764B", "burlywood", "tan3")

relDendA = ggplot() +
  geom_rect(aes(xmin = 21.25, xmax = 52.75 , ymin = 0.475, ymax = 0.69), fill = "black", alpha = 0.15) +
  geom_rect(aes(xmin = 53.25, xmax = 70.75 , ymin = 0.425, ymax = 0.75), fill = "black", alpha = 0.15) +
  geom_segment(data = segment(relateDData), aes(x = x, y = 1-y, xend = xend, yend = 1-yend2), size = 0.5) +
  geom_point(data = relDendPoints2, aes(x = x, y = 1-y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
  geom_text(data = relDendPoints2, aes(x = x, y = (1 - y + .05), label = label), angle = 90) +
  geom_text(aes(x = 37, y = 0.45), label = "Sib 1", size = 6) + 
  geom_text(aes(x = 62, y = 0.40), label = "Sib 2", size = 6) + 
  #scale_fill_brewer(palette = "Dark2", name = "Population") +
  scale_fill_manual(values = flPal, name= "Site:")+
  scale_shape_manual(values = c(24, 25), name = "Depth Zone:")+
 # spacing technical replicates further from leaf
 # geom_text(data = dendPoints, aes(x = x, y = (y - .025), label = label), angle = 90) +
  labs(y = "Relatedness") +
  scale_y_reverse() +
  guides(fill = guide_legend(override.aes = list(shape = 22, size = 10), ncol = 2), shape = guide_legend(override.aes = list(size = 8), ncol = 1, order = 1)) +
  coord_cartesian(xlim = c(5, 215)) +
  theme_classic()

relDend = relDendA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 24, color = "black", angle = 90),
  axis.text.y = element_text(size = 20, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 24),
  legend.text = element_text(size = 20),
  legend.position = "bottom")

# relDend

Now plotting the same dendrogram colored by Admixture group

relDendKA = ggplot() +
  geom_rect(aes(xmin = 21.25, xmax = 52.75 , ymin = 0.475, ymax = 0.69), fill = "black", alpha = 0.15) +
    geom_rect(aes(xmin = 53.25, xmax = 70.75 , ymin = 0.425, ymax = 0.75), fill = "black", alpha = 0.15) +
  geom_segment(data = segment(relateDData), aes(x = x, y = 1-y, xend = xend, yend = 1-yend2), size = 0.5) +
  geom_point(data = relDendPoints2, aes(x = x, y = 1-y, fill = cluster, shape = depth), size = 4, stroke = 0.25) +
  geom_text(data = relDendPoints2, aes(x = x, y = (1 - y + .05), label = label), angle = 90) +
  geom_text(aes(x = 37, y = 0.45), label = "Sib 1", size = 6) + 
  geom_text(aes(x = 62, y = 0.40), label = "Sib 2", size = 6) + 
  #scale_fill_brewer(palette = "Dark2", name = "Population") +
  scale_fill_manual(values = c(kColPal, "gray50"), name= "Lineage:") +
  scale_shape_manual(values = c(24, 25), name = "Depth Zone:") +
 # spacing technical replicates further from leaf
 # geom_text(data = dendPoints, aes(x = x, y = (y - .025), label = label), angle = 90) +
  labs(y = "Relatedness") +
  scale_y_reverse() +
  guides(fill = guide_legend(override.aes = list(shape = 22, size = 10), ncol = 3), shape = guide_legend(override.aes = list(size = 8), ncol = 1, order = 1)) +
  coord_cartesian(xlim = c(5,215)) +
  theme_classic()

relDendK = relDendKA + theme(
  axis.title.x = element_blank(),
  axis.text.x = element_blank(),
  axis.line.x = element_blank(),
  axis.ticks.x = element_blank(),
  axis.title.y = element_text(size = 24, color = "black", angle = 90),
  axis.text.y = element_text(size = 20, color = "black"),
  axis.line.y = element_line(),
  axis.ticks.y = element_line(),
  panel.grid = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  legend.key = element_blank(),
  legend.title = element_text(size = 24),
  legend.text = element_text(size = 20),
  legend.position = "bottom")

# relDendK

Put both plots together with Patchwork

relDendPlots = (relDend / relDendK) +
  plot_annotation(tag_levels = 'A') +
  plot_layout(guides = "collect") & 
  theme(plot.tag = element_text(size = 32),
        legend.position = "bottom")

ggsave("../figures/figureS2.png", plot = relDendPlots, height = 13, width = 35, units = "in", dpi = 300)
ggsave("../figures/figureS2.svg", plot = relDendPlots, height = 13, width = 35, units = "in", dpi = 300)


Heterozygosity, inbreeding, and relatedness plots

Heterozygosity across all RAD loci

hetMelt = melt(het, id.vars = c("sample", "Site", "Depth"), variable.name = "type", value.name = "heterozygosity")

hetMelt = melt(domClust, id.vars = c("sample", "Site", "Depth"), variable.name = "type", value.name = "heterozygosity")

hetTheme = theme(axis.title.x = element_blank(),
        axis.text.x = element_text(color = "black", size = 12),
        axis.ticks.x = element_line(color = "black"),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 12),
        axis.ticks.y = element_line(color = "black"),
        legend.key.size = unit(0.3, 'cm'),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_rect(color = "black"),
        panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

dodge <- position_dodge(width = 1)

hetPlotAll = ggplot(data = het,aes(x = Depth, y = All, fill = Site)) +
  geom_point(aes(color = Site),shape = 16, position = position_jitterdodge(seed = 1, dodge.width = 1), size = 0.75, alpha = 0.4) +
  geom_violin(size = 0.4, alpha = 0.5, position = dodge, color = "gray45") +
  annotate(geom = "text", x = 2.20, y = 0.009, label = bquote(~bold("Depth:")), size = 3) +
  annotate(geom = "text", x = 2.3, y = 0.0087, label = bquote(~italic(F)~" = "~.(round(hetAllANOVA[[1]][2,4],3))), size = 3) +
  annotate(geom = "text", x = 2.26, y = 0.0084, label = bquote(~italic(p)~" < 0.01"), size = 3) +
  guides(color = guide_legend(override.aes = list(shape = 22, size = 3, fill = flPal)), fill = "none", shape = "none", size = "none") +
  xlab("Population depth zone") +
  ylab("Heterozygosity (All RAD loci)") +
  scale_fill_discrete(type = rev(flPal), breaks=c('Upper Keys', 'Lower Keys', 'Tortugas Bank', "Riley's Hump")) +  
  scale_color_discrete(type = rev(flPal), breaks=c('Upper Keys', 'Lower Keys', 'Tortugas Bank', "Riley's Hump")) +
  theme_bw() +
  theme(legend.position = c(0.25, 0.85),
        legend.background = element_blank()) +
  hetTheme 

Heterozygosity across SNP loci

hetPlotSnps = ggplot(data = het, aes(x = Depth, y = SNPs, fill = Site)) +
  geom_point(aes(color = Site),shape = 16, position = position_jitterdodge(seed = 1, dodge.width = 1), size = 0.75, alpha = 0.4) +
  geom_violin(size = 0.4, alpha = 0.5, position = dodge, color = "gray45") +
  annotate(geom = "text", x = 2.16, y = 0.475, label = bquote(~bold("Depth:")), size = 3) +
    annotate(geom = "text", x = 2.28, y = 0.455, label = bquote(~italic(F)~" = "~.(round(hetSnpANOVA[[1]][2,4],3))), size = 3) +
  annotate(geom = "text", x = 2.26, y = 0.435, label = bquote(~italic(p)~" < 0.0001"), size = 3) +
  xlab("Population depth zone") +
  ylab("Heterozygosity (SNP loci)") +
  scale_fill_manual(values = rev(flPal)) +
  scale_y_continuous(breaks=seq(0.1, 0.5, by = .05)) +
  scale_color_manual(values = rev(flPal)) +
  theme_bw() +
  theme(legend.position = "none") +
  hetTheme

Mean inbreeding plot

inbreedingPlot = ggplot(data = het, aes(x = Depth, y = inbreed, fill = Site)) +
  geom_point(aes(color = Site),shape = 16, position = position_jitterdodge(seed = 1, dodge.width = 1), size = 0.75, alpha = 0.4) +
  geom_violin(size = 0.4, alpha = 0.5, position = dodge, color = "gray45") +
  annotate(geom = "text", x = 2.16, y = 0.35, label = bquote(~bold("Depth:")), size = 3) +
   annotate(geom = "text", x = 2.28, y = 0.33, label = bquote(~italic(F)~" = "~.(round(inbreedANOVA[[1]][2,4],3))), size = 3) +
  annotate(geom = "text", x = 2.26, y = 0.31, label = bquote(~italic(p)~" < 0.0001"), size = 3) +
  xlab("Population depth zone") +
  ylab(bquote(~"Inbreeding coefficient ("*italic(F)*")")) +
  scale_fill_manual(values = rev(flPal)) +  
  scale_color_manual(values = rev(flPal)) +
  scale_y_continuous(breaks=seq(0, 0.4, by = .05)) +
  theme_bw() +
  theme(legend.position = "none") +
  hetTheme

# inbreedingPlot

Pairwise relatedness plot

relatePlot = ggplot(data = sintRelate, aes(x = Depth, y = Rab, fill = Site)) +
  geom_point(aes(color = Site),shape = 16, position = position_jitterdodge(seed = 1, dodge.width = 1), size = 0.75, alpha = 0.4) +
  geom_violin(size = 0.25, alpha = 0.5, position = dodge, color = "gray35") +
  annotate(geom = "text", x = 2.16, y = 0.51, label = bquote(~bold("Depth:")), size = 3) +
  annotate(geom = "text", x = 2.28, y = 0.48, label = bquote(~italic(F)~" = "~.(round(relateANOVA[[1]][2,4],3))), size = 3) +
  annotate(geom = "text", x = 2.26, y = 0.45, label = bquote(~italic(p)~" < 0.0001"), size = 3) +
  xlab("Population depth zone") +
  ylab(bquote(~"Relatedness ("*italic(R)["AB"]*")")) +
  scale_fill_manual(values = rev(flPal)) +
  scale_color_manual(values = rev(flPal)) +
  scale_y_continuous(breaks=seq(0, 0.5, by = .1)) +
  guides() +
  theme_bw() +
  theme(legend.position = "none") +
  hetTheme

Finally, put all plots together with Patchwork

hetPlots = ((hetPlotAll | hetPlotSnps) / (inbreedingPlot | relatePlot)) +
  plot_annotation(tag_levels = 'A')&
  theme(plot.tag = element_text(size = 16))

# hetPlots

ggsave("../figures/figure2.png", plot = hetPlots, width = 18, height = 17, units = "cm", dpi = 300)
ggsave("../figures/figure2.svg", plot = hetPlots, width = 18, height = 17, units = "cm", dpi = 300)

Outliers

bayescan = read.table("../data/snps/fkSint.baye_fst.txt",header=T) %>% mutate(loc = rownames(.), out.05 = ifelse(qval < 0.05, 1, 0), out.1 = ifelse(qval < 0.1, 1, 0))
bayescan[bayescan[, 3]<=0.0001, 3] = 0.0001

bayeScEnv = read.table("../data/snps/fkSint.bayeS_fst.txt", header=T) %>% filter(qval_g < 0.05) %>% mutate(loc = rownames(.), depthOut = 1) %>% dplyr::select(loc, depthOut)

bayescan = bayescan %>% left_join(bayeScEnv) 
## Joining, by = "loc"
bayescan$depthOut = bayescan$depthOut %>% replace_na(0)

sum(bayescan$out.05)
## [1] 591
sum(bayescan$out.1)
## [1] 762
sum(bayescan$depthOut)
## [1] 6
for(i in 1:nrow(bayescan)){
  if(bayescan$depthOut[i] == 1){
    bayescan$out.05[i] = 2
  }
}

bayescanPlotA = ggplot(data = bayescan, aes(x = log10(qval), y = fst, color = as.factor(out.05), alpha = as.factor(out.05))) +  
  geom_point(size = 1) +
  geom_vline(xintercept = log10(0.05), linetype = 2, color = purple) +
  xlab(expression(log[10]*"("*italic("q")*"-value)")) +
  ylab(expression(italic("F")[ST])) +
  scale_x_reverse() +
  scale_color_manual(values = c("grey45", purple, pink)) +
  scale_alpha_manual(values = c(0.25, 0.25, 0.5)) +
  theme_bw()

bayescanPlot = bayescanPlotA +
        theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_text(color = "black", size = 12),
        axis.ticks.x = element_line(color = "black"),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        axis.ticks.y = element_line(color = "black"),
        legend.position = "none",
        legend.key.size = unit(0.3, 'cm'),
        panel.border = element_rect(color = "black"),
        panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# bayescanPlot

ggsave("../figures/figureS3.svg", plot = bayescanPlot, width = 14, height = 8, units = "cm", dpi = 300)
ggsave("../figures/figureS3.png", plot = bayescanPlot, width = 14, height = 8, units = "cm", dpi = 300)

Now we can look at the minor allele frequencies of the outliers by population

fkSintAF = read.delim("../data/snps/fkSintAlleleFreq.txt")

fkSintAF$chromLoc = paste(fkSintAF$chrom,fkSintAF$pos)

fkSintAF$pop = factor(fkSintAF$pop)
fkSintAF$pop = factor(fkSintAF$pop, levels = levels(fkSintAF$pop)[c(3, 4, 5, 6, 1, 2, 7, 8)], labels = c("Riley's Hump\nMesophotic", "Riley's Hump\nShallow", "Tortugas Bank\nMesophotic", "Tortugas Bank\nShallow", "Lower Keys\nMesophotic", "Lower Keys\nShallow", "Upper Keys\nMesophotic", "Upper Keys\nShallow"))

fkSintAF = fknmsSites[-c(66,68,164,166,209,211),] %>%  group_by(site, depthZone) %>% summarise(avgDepthM = mean(depthM)) %>% mutate(pop = paste(site, depthZone, sep = "\n")) %>% dplyr::select(pop, site, depthZone, avgDepthM) %>% left_join(fkSintAF, .)
## `summarise()` has grouped output by 'site'. You can override using the `.groups`
## argument.
## Joining, by = "pop"
fkSintAF$pop = factor(fkSintAF$pop)
fkSintAF$pop = factor(fkSintAF$pop, levels = levels(fkSintAF$pop)[c(3, 4, 5, 6, 1, 2, 7, 8)])

mafTheme = theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        axis.ticks.y = element_blank(),
        legend.position = "right",
        legend.key.size = unit(0.3, 'cm'),
        legend.title = element_text(color = "black", size = 12),
        legend.text = element_text(color = "black", size = 10),
        panel.border = element_rect(color = "black", size = .25),
        panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

mafPlot = ggplot(data = fkSintAF, aes(x = chromLoc, y = pop, fill = minFreq)) +
  geom_tile() +
  geom_segment(data = fkSintAF, aes(x = -7, xend = -220, y = pop, yend = pop, color = site), size = 8.5) +
  scale_color_manual(values = flPal, guide = NULL) +
  scale_fill_gradientn(colours = paletteer_d("RColorBrewer::RdPu"), limit = c(0, 1), name = "Minor allele frequency")+
  xlab("Outlier locus") +
  ylab("Population") +
  scale_y_discrete() +
  coord_cartesian(xlim = c(-1, 592), clip = "off") +
  guides(fill = guide_colorbar(title.position = "left", barheight = unit(100, "pt"))) + 
  theme_bw() +
  mafTheme +
  theme(axis.text.y = element_text(color = "gray90", hjust = 1, size = 9),
        legend.title = element_text(angle = 90))

# mafPlot

Linear regression of minor allele frequencies by depth

mafLm = summary(lm(data = fkSintAF, minFreq ~ avgDepthM))
r2 = round(mafLm$adj.r.squared, 3)

fkSintAF$depthZone = factor(fkSintAF$depthZone, levels = levels(fkSintAF$depthZone)[c(2, 1)])

mafLrPlot = ggplot(data = fkSintAF, aes(x = avgDepthM, y = minFreq)) +
  # geom_point(aes(color = site, fill = site),shape = 21, position = position_jitter(width = 0.5), size = 0.5, alpha = 0.4) +
  geom_violin(aes(group = pop, fill = site, linetype = depthZone, color = site), size = 0.25, width = 2, alpha = 0.5, scale = "width") +
  # geom_boxplot(aes(group = pop, fill = site, linetype = depthZone), size = 0.25, alpha = 0.5, position = dodge, width = 1.5,color = "gray35", outlier.colour = NA) +
  geom_smooth(method = "lm", color = pink, size = 0.25) +
  annotate(geom = "text", x = 37, y = -0.045, label = bquote(~italic(R^2)~" = "~.(r2)*", "~italic(p)~" < 0.001"), size = 3) +
  xlab("Population depth (m)") +
  ylab("Outlier locus \nminor allele frequency") +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = flPal, name = "Site") +
  scale_linetype_discrete(name = "Depth zone", breaks = c("Shallow", "Mesophotic")) +
  theme_bw() +
  mafTheme +
  theme(axis.text.x = element_text(color = "black", size = 10),
        axis.ticks.x = element_line(color = "black"))

# mafLrPlot

Inbred outliers

Subsetting the inbred shallow individuals to look at their minor allele frequencies vs other shallow individuals

shalInbred = sintBreed
shalInbred$a = c(1:220)
bams = read.delim("../data/snps/bamsNoClones", header = FALSE) %>% mutate(a = c(1:220))
pops = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("Site" = site, "Depth" = depthZone) %>% mutate(a = c(1:220), pop = paste(Site, Depth, sep = " ")) %>% mutate(across(where(is.character), str_remove_all, pattern = fixed("'")))  %>% 
  mutate(across(where(is.character), str_remove_all, pattern = fixed(" ")))# Reads in population data

shalBams = bams %>% left_join(pops) %>% left_join(shalInbred) %>% dplyr::select(-a) %>% filter(Depth == "Shallow")
## Joining, by = "a"
## Joining, by = "a"
shalBams$pop = factor(shalBams$pop)

head(shalBams)
##                        V1         Site   Depth                 pop         inbreed
## 1 fk_S010.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.1193433952381
## 2 fk_S012.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.0409269615385
## 3 fk_S013.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.0561628115942
## 4 fk_S016.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.0000005147059
## 5 fk_S026.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.0712638659794
## 6 fk_S027.trim.un.bt2.bam TortugasBank Shallow TortugasBankShallow 0.3228964818653
shalInbredBams =  shalBams %>% filter(inbreed >= 0.25) %>% dplyr::select(V1, pop) %>% group_split(pop) %>% setNames(unique(shalBams$pop))

shalOutbredBams = shalBams %>% filter(inbreed < 0.25) %>% dplyr::select(V1, pop) %>% group_split(pop) %>% setNames(unique(shalBams$pop))

write.table(shalOutbredBams$UpperKeysShallow, "../data/snps/shalUKOutbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalOutbredBams$LowerKeysShallow, "../data/snps/shalLKOutbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalOutbredBams$TortugasBankShallow, "../data/snps/shalTBOutbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalOutbredBams$RileysHumpShallow, "../data/snps/shalRHOutbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)

write.table(shalInbredBams$UpperKeysShallow, "../data/snps/shalUKInbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalInbredBams$LowerKeysShallow, "../data/snps/shalLKInbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalInbredBams$TortugasBankShallow, "../data/snps/shalTBInbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table(shalInbredBams$RileysHumpShallow, "../data/snps/shalRHInbred.pop", row.names = FALSE, col.names = FALSE, quote = FALSE)

Use these files to split outlier .vcf into inbred//outbred groups in each shallow population

shalInbred = shalBams
shalInbred$inbred = NA

for(i in 1:nrow(shalInbred)) {
  if(shalInbred$inbreed[i] >= 0.25) {
   shalInbred$inbred[i] = "Inbred"
  } else {
    shalInbred$inbred[i] = "Outbred"
  }
}

shalInbred = shalInbred %>% mutate(pop = paste(Site, inbred)) 
shalInbred = shalInbred %>% group_by(pop) %>% summarise (avgInbreed = mean(inbreed)) 
shalInbred$pop = c("Lower Keys Inbred", "Lower Keys Outbred", "Riley's Hump Inbred", "Riley's Hump Outbred", "Tortugas Bank Inbred", "Tortugas Bank Outbred", "Upper Keys Inbred", "Upper Keys Outbred")

fkSintInbredAF = read.delim("../data/snps/fkSintInbredAlleleFreq.txt")

fkSintInbredAF$chromLoc = paste(fkSintInbredAF$chrom,fkSintInbredAF$pos)

fkSintInbredAF$pop = factor(fkSintInbredAF$pop)
fkSintInbredAF$pop = factor(fkSintInbredAF$pop, levels = levels(fkSintInbredAF$pop)[c(3, 4, 5, 6, 1, 2, 7, 8)], labels = c("Riley's Hump", "Riley's Hump", "Tortugas Bank", "Tortugas Bank", "Lower Keys", "Lower Keys", "Upper Keys", "Upper Keys"))

fkSintInbredAF = fkSintInbredAF %>% mutate(site = pop) %>% dplyr::select(-pop)
fkSintInbredAF = fkSintInbredAF %>% mutate(pop = paste(site, inbred)) 
fkSintInbredAF = fkSintInbredAF %>% left_join(shalInbred)
## Joining, by = "pop"
fkSintInbredAF$pop = factor(fkSintInbredAF$pop)
fkSintInbredAF$pop = factor(fkSintInbredAF$pop, levels = levels(fkSintInbredAF$pop)[c(4, 3, 6, 5, 2, 1, 8, 7)], labels = c("Riley's Hump\nOutbred", "Riley's Hump\nInbred", "Tortugas Bank\nOutbred", "Tortugas Bank\nInbred", "Lower Keys\nOutbred", "Lower Keys\nInbred", "Upper Keys\nOutbred", "Upper Keys\nInbred"))

fkSintInbredAF$site = factor(fkSintInbredAF$site)
fkSintInbredAF$site = factor(fkSintInbredAF$site, levels = levels(fkSintInbredAF$site)[c(4,3,2,1)])

fkSintInbredAF$inbred = factor(fkSintInbredAF$inbred)
fkSintInbredAF$inbred = factor(fkSintInbredAF$inbred, levels = levels(fkSintInbredAF$inbred)[c(2,1)])

mafInbredPlot = ggplot(data = fkSintInbredAF, aes(x = chromLoc, y = pop, fill = minFreq)) +  
  geom_tile() +
  geom_segment(data = fkSintInbredAF, aes(x = -7, xend = -220, y = pop, yend = pop, color = site), size = 8.5) +
  scale_fill_gradientn(colours = paletteer_d("RColorBrewer::RdPu"), limit = c(0, 1), name = "Minor allele frequency")+
  scale_color_manual(values = flPal, guide = NULL) +
  xlab("Outlier locus") +
  ylab("Shallow population") +
  coord_cartesian(xlim = c(-1, 592), clip = "off", expand = TRUE) +
  guides(fill = guide_colorbar(title.position = "left", barheight = unit(100, "pt"))) +
  theme_bw() +
  mafTheme +
  theme(axis.text.y = element_text(color = "gray90", hjust = 1, size = 9),
        legend.title = element_text(angle = 90))

# mafInbredPlot
mafInbredLm = summary(lm(data = fkSintInbredAF, minFreq ~ avgInbreed))
inbredR2 = round(mafInbredLm$adj.r.squared, 3)

mafLrInbredPlot = ggplot(data = fkSintInbredAF, aes(x = avgInbreed, y = minFreq)) +
  # geom_point(aes(color = site, fill = site),shape = 21, position = position_jitter(width = 0.005), size = 0.5, alpha = 0.4) +
  geom_violin(aes(group = pop, fill = site, linetype = inbred, color = site), size = 0.25, alpha = 0.5, width = 0.025, scale = "width") +
  # geom_boxplot(aes(group = pop, fill = site, linetype = inbred), size = 0.25, alpha = 0.5, position = dodge, width = .015, color = "gray35", outlier.colour = NA) +
  geom_smooth(method = "lm", color = pink, size = 0.25) +
  annotate(geom = "text", x = 0.26, y = -0.04, label = bquote(~italic(R^2)~" = "~.(inbredR2)*", "~italic(p)~" < 0.001"), size = 3) +
  xlab("Mean inbreeding coefficient") +
  ylab("Outlier locus \nminor allele frequency") +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = flPal, name = "Site") +
  scale_linetype_discrete(name = "Inbreeding", breaks = c("Inbred", "Outbred")) +
  theme_bw() +
  mafTheme +
  theme(axis.text.x = element_text(color = "black", size = 10),
        axis.ticks.x = element_line(color = "black"))

# mafLrInbredPlot
mafPlots = (mafPlot | mafLrPlot) / (mafInbredPlot | mafLrInbredPlot) +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 16))

ggsave("../figures/figure3.png", plot = mafPlots, width = 25, height = 15, units = "cm", dpi = 300)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave("../figures/figure3.svg", plot = mafPlots, width = 25, height = 15, units = "cm", dpi = 300)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'


Population structure

Discriminant analysis of principal components

dapcKPal = paletteer_d("wesanderson::Darjeeling1")[c(1,4,2,5)]
flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]

sintVcf = read.vcfR("../data/snps/sintNoClones.bcf", verbose = FALSE)
sintGenlightPopDepth = vcfR2genlight(sintVcf, n.cores = 2) # Converts the vcf file into a file format that poppr uses the "genlight" format

locNames(sintGenlightPopDepth) = paste(sintVcf@fix[,1],sintVcf@fix[,2],sep="_")
popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone, "depthM" = depthM) # Reads in population data for each sample
popData$popdepth = paste(popData$pop, popData$depth, sep = " ")

popData$depth = factor(popData$depth)
popData$depth = factor(popData$depth, levels(popData$depth)[c(2,1)])

popData$pop = factor(popData$pop)
popData$pop = factor(popData$pop, levels(popData$pop)[c(4, 1, 3, 2)])

popData$popdepth = factor(popData$popdepth)
popData$popdepth = factor(popData$popdepth, levels(popData$popdepth)[c(8,7,2,1,6,5,4,3)])
levels(popData$popdepth) = c("UK Shallow", "UK Mesophotic",  "LK Shallow",  "LK Mesophotic", "TB Shallow", "TB Mesophotic", "RH Shallow", "RH Mesophotic")

strata(sintGenlightPopDepth) = data.frame(popData)
setPop(sintGenlightPopDepth) = ~popdepth

sintGenlightPopDepth$pop = factor(sintGenlightPopDepth$pop)
sintGenlightPopDepth$pop = factor(sintGenlightPopDepth$pop, 
                                levels(sintGenlightPopDepth$pop)[c(7, 8, 6, 5, 2, 1, 4, 3)])

Find the best number of clusters. You may need to run WITHOUT [choose.n.clust = FALSE] and [criterion = …] In this case you would examine the BIC plot and pick nClust interactively in the console

grp = find.clusters(sintGenlightPopDepth, n.pca = 219, max.n.clust = 11, choose.n.clust = FALSE, criterion = "goesup")

#4 clusters found
kVals = data.frame(K = c(1:11),BIC = grp$Kstat)

bicPlotA = ggplot(data = kVals, aes (x = K, y = BIC)) +
                 geom_line(color = "gray65") +
                 geom_point(color = "gray65", size = 1) +
                 geom_point(x = 4, y = grp$stat, color = pink, size = 1) +
                 xlab(expression(paste("Number of clusters (",italic(K),")",sep = ""))) +
                 ylab("BIC") +
                 coord_cartesian(ylim = c(1610,1640)) +
                 scale_x_continuous(breaks = c(seq(1,11,2))) +
                 scale_y_continuous(breaks = c(seq(1610,1640,10))) +
                 theme_bw()

bicPlot = bicPlotA +
        theme(axis.title.x = element_text(color = "black", size = 8),
        axis.text.x = element_text(color = "black", size = 6),
        axis.ticks.x = element_line(color = "black"),
        axis.title.y = element_text(color = "black", size = 8),
        axis.text.y = element_text(color = "black", size = 6),
        axis.ticks.y = element_line(color = "black"),
        legend.position = "none",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_blank(),
        plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

bicPlot

table(pop(sintGenlightPopDepth), grp$grp)
##                
##                  1  2  3  4
##   UK Shallow    15  5  0 10
##   UK Mesophotic 27  3  0  0
##   LK Shallow    11  7  3  9
##   LK Mesophotic 29  1  0  0
##   TB Shallow     9  8 10  3
##   TB Mesophotic 19  6  0  0
##   RH Shallow    16  5  2  7
##   RH Mesophotic  8  5  0  2
table.value(table(pop(sintGenlightPopDepth), grp$grp), col.lab=paste("inf", 1:4))

set.seed(694)
xval1A = xvalDapc(tab(sintGenlightPopDepth, NA.method = "mean"), grp$grp, n.rep = 100)

xval1A$`Number of PCs Achieving Lowest MSE`
## [1] "60"
# Running more replicates around the peak of the last cross validation (~80)
# This takes ~6 hrs (!!!)
# I have saved the results in .RData to continue using, but if you need to run with new data, un-comment and run the following

# set.seed(694)
# system.time(xval1 <- xvalDapc(tab(sintGenlightPopDepth, NA.method = "mean"), grp$grp, n.pca = 50:120, n.rep = 1000, parallel = "multicore", ncpus = 4))

xval1$`Number of PCs Achieving Highest Mean Success`
## [1] "55"
# 55 PCs retained
 
dapc1 = dapc(sintGenlightPopDepth, grp$grp, n.pca = as.numeric(xval1$`Number of PCs Achieving Highest Mean Success`), n.da = 4)

dapc1DF = as.data.frame(dapc1$ind.coord) %>% cbind(grp = dapc1$grp)

scatter(dapc1)

df1PlotA = ggplot(data = dapc1DF, aes(x = LD1, group = grp, fill = grp, color = grp)) + 
  geom_point(y = 0, shape = "|", size = 3, alpha = 0.7) + 
  geom_density(alpha = 0.7, size = 0.25) +
  xlab("Discriminant function 1") +
  ylab("Density") +
  scale_color_manual(values = dapcKPal, name = "Cluster") +
  scale_fill_manual(values = dapcKPal, name = "Cluster") +
  guides(color = "none", fill = "none", shape = "none") +
  theme_bw()

df1Plot = df1PlotA +
        theme(axis.title.x = element_text(color = "black", size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 12),
        axis.text.y = element_text(color = "black", size = 10),
        axis.ticks.y = element_line(color = "black"),
        legend.position = "bottom",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_blank(),
        plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

df1Plot

set.seed(694)
xval2A <- xvalDapc(tab(sintGenlightPopDepth, NA.method = "mean"), pop(sintGenlightPopDepth), n.rep = 100)

xval2A$`Number of PCs Achieving Lowest MSE`
## [1] "60"
# 60

# Running more replicates around the peak of the last cross validation (60)
# This takes a long time! I'm commenting it out and loading it in with the .RData file for this project

# set.seed(694)
# xval2 <- xvalDapc(tab(sintGenlightPopDepth, NA.method = "mean"), pop(sintGenlightPopDepth), n.pca = 30:90, n.rep = 1000, parallel = "multicore", ncpus = 3)

xval2$`Number of PCs Achieving Lowest MSE`
## [1] "61"
# 61 PCs retained

dapc2 = dapc(sintGenlightPopDepth, sintGenlightPopDepth$pop, n.pca = as.numeric(xval2$`Number of PCs Achieving Lowest MSE`), n.da = 7)

sintDapcVar = round(dapc2$eig/sum(dapc2$eig)*100, 1)
head(sintDapcVar)
## [1] 28.7 23.9 13.6 12.7 10.5  6.6
grpCentroids = as.data.frame(dapc2$grp.coord)
colnames(grpCentroids) = c("popLD1", "popLD2", "popLD3", "popLD4", "popLD5", "popLD6", "popLD7")
grpCentroids$popdepth = levels(popData$popdepth)

sintDapc = popData %>% cbind(dapc2$ind.coord) %>% left_join(grpCentroids) %>% cbind(K = dapc1$grp)
## Joining, by = "popdepth"
sintDapcPlotA = ggplot(sintDapc, aes(x = LD1, y = LD2, fill = pop, shape = depth)) +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  geom_segment(aes(x = popLD1, y = popLD2, xend = LD1, yend = LD2, color = pop), alpha = 0.5) +
  geom_point(aes(x = LD1, y = LD2, shape = depth, color = pop), size = 2, alpha = 0.5, show.legend = TRUE) +
  scale_shape_manual(values = c(24, 25), name = "Depth Zone") +
  # stat_ellipse(data = sintDapc, aes(color = pop, linetype = depth), type = "t"  , geom = "polygon", alpha = 0.1, size = 0.4) +
  geom_point(aes(x = popLD1, y = popLD2, shape = depth), size = 4, color = "gray15") +
  # geom_label(aes(x = popLD1, y = popLD2, label = popdepth, fill = pop), size = 2, color = "gray90") +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = flPal, guide = NULL) +
  xlab(paste("DA 1 (", sintDapcVar[1], "%)", sep = "")) + 
  ylab(paste("DA 2 (", sintDapcVar[2], "%)", sep = "")) +
  # guides(shape = "none", fill = guide_legend(override.aes = list(shape = 22, color = NA, alpha = NA), order = 1)) +
  guides(shape = guide_legend(override.aes = list(size = 3)), fill = guide_legend(override.aes = list(shape = 22, color = NA, alpha = NA), order = 1)) +
  theme_bw()

sintDapcPlot = sintDapcPlotA +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.position = "right",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

sintDapcPlot

# ggsave("../figures/dapc.png", plot = sintDapcPlot, height = 12, width = 18, units = "cm", dpi = 300)
sintDapcKPlotA = ggplot(sintDapc, aes(x = LD1, y = LD2, fill = K, shape = depth)) +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  geom_point(aes(x = LD1, y = LD2, shape = depth, color = K), size = 2, alpha = 0.5) +
  scale_shape_manual(values = c(24,25), name = "Depth Zone") +
  scale_fill_manual(values = dapcKPal, name = "Cluster") +
  scale_color_manual(values = dapcKPal, name = "Cluster", guide = NULL) +
  xlab(paste("DA 1 (", sintDapcVar[1], "%)", sep = "")) + 
  ylab(paste("DA 2 (", sintDapcVar[2], "%)", sep = "")) +
  guides(shape = guide_legend(override.aes = list(size = 3)), fill = guide_legend(override.aes = list(shape = 22,  alpha = NA, color = NA, size = 4), order = 1))+
  theme_bw()

sintDapcKPlot = sintDapcKPlotA +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.position = "right",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

sintDapcKPlot

sintDapcDepthPlotA = ggplot(sintDapc, aes(x = LD1, y = LD2, fill = depthM, shape = depth)) +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  # geom_segment(aes(x = popLD1, y = popLD2, xend = LD1, yend = LD2, color = pop)) +
  geom_point(aes(x = LD1, y = LD2, shape = depth), color = "black", size = 2.5, alpha = 0.7, show.legend = TRUE) +
  scale_shape_manual(values = c(24,25), name = "Depth Zone") +
  # stat_ellipse(data = sintDapc, type = "t"  , geom = "polygon", alpha = 0.1, fill = NA, size = 0.4) +
  # geom_point(aes(x = popLD1, y = popLD2, shape = depth), size = 5, color = "black") +
  # geom_label(aes(x = popLD1, y = popLD2, label = popdepth), fill = "white", size = 2) +
  # scale_fill_gradient(limit = c(10, 50), low = "#39E6F4", high = "#002C7C", name = "Depth (m)") +
  scale_fill_gradientn(colours = brewer.pal(11, "RdBu"), limit = c(14, 46), name = "Depth (m)", breaks = c(15, 20, 25, 30, 35, 40, 45)) +
  xlab(paste("DA 1 (", sintDapcVar[1], "%)", sep = "")) + 
  ylab(paste("DA 2 (", sintDapcVar[2], "%)", sep = "")) +
  guides(shape = guide_legend(override.aes = list(size = 3.5, color = "black", fill = "white", stroke = 0.5, alpha = 0.7), order = 2), fill = guide_colorbar(reverse = TRUE, order = 1, barwidth = 0.7)) +
  theme_bw()

sintDapcDepthPlot = sintDapcDepthPlotA +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.position = "right",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

sintDapcDepthPlot


Analysis of Molecular Variance

sintVcf = read.vcfR("../data/snps/sintNoClones.bcf", verbose = FALSE)

sintGenlightPopDepth = vcfR2genlight(sintVcf, n.cores = 2) # Converts the vcf file into a file format that poppr uses the "genlight" format

locNames(sintGenlightPopDepth) = paste(sintVcf@fix[,1],sintVcf@fix[,2],sep="_")
popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone) # Reads in population data for each sample

popData$popdepth = paste(popData$pop, popData$depth, sep = " ")

strata(sintGenlightPopDepth) = data.frame(popData)
setPop(sintGenlightPopDepth) = ~popdepth

#Run Analysis of MOlecular VAriance (AMOVA)
amova <- poppr.amova(sintGenlightPopDepth, ~popdepth) 
amova
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                                  Df     Sum Sq  Mean Sq
## Between popdepth                  7   35867.77 5123.967
## Between samples Within popdepth 212  680293.73 3208.933
## Within samples                  220  562159.50 2555.270
## Total                           439 1278321.00 2911.893
## 
## $componentsofcovariance
##                                                  Sigma         %
## Variations  Between popdepth                  34.98402   1.19928
## Variations  Between samples Within popdepth  326.83112  11.20403
## Variations  Within samples                  2555.27045  87.59669
## Total variations                            2917.08559 100.00000
## 
## $statphi
##                            Phi
## Phi-samples-total    0.1240331
## Phi-samples-popdepth 0.1134003
## Phi-popdepth-total   0.0119928
#Calculate significance levels of the AMOVA with 999 permutations; this takes a while 
set.seed(694)
amovasignif <- randtest(amova, nrepet = 999)

amovaPerc = paste(round(amova$componentsofcovariance$`%`[1], 2), "%",sep="")
amovaP = amovasignif$pvalue[3]

paste(amovaPerc, amovaP)
## [1] "1.2% 0.001"


Population differentiation

Calculating differentiation with Pairwise Fst

sintVcf = read.vcfR("../data/snps/sintNoClones.bcf", verbose = FALSE)
sintGenlightPopDepth = vcfR2genlight(sintVcf, n.cores = 2) # Converts the vcf file into a file format that poppr uses the "genlight" format
locNames(sintGenlightPopDepth) = paste(sintVcf@fix[,1],sintVcf@fix[,2],sep="_")
popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone) # Reads in population data for each sample
popData$popdepth = paste(popData$pop, popData$depth, sep = " ")

strata(sintGenlightPopDepth) = data.frame(popData)
setPop(sintGenlightPopDepth) = ~popdepth

sintGenlightPopDepth$pop = factor(sintGenlightPopDepth$pop)
sintGenlightPopDepth$pop = factor(sintGenlightPopDepth$pop, 
                                levels(sintGenlightPopDepth$pop)[c(7, 8, 6, 5, 2, 1, 4, 3)])
                                
levels(sintGenlightPopDepth$pop) = c("Upper Keys\nShallow", "Upper Keys\nMesophotic", "Lower Keys\nShallow", "Lower Keys\nMesophotic", "Tortugas Bank\nShallow", "Tortugas Bank\nMesophotic", "Riley's Hump\nShallow", "Riley's Hump\nMesophotic")

set.seed(694)
fk.fst <- stamppFst(sintGenlightPopDepth, nboots = 99, percent = 95, nclusters = 2) 
fk.fst$Fsts
##                           Tortugas Bank\nMesophotic Tortugas Bank\nShallow
## Tortugas Bank\nMesophotic                        NA                     NA
## Tortugas Bank\nShallow                  0.021928230                     NA
## Riley's Hump\nMesophotic                0.002397382            0.012550165
## Riley's Hump\nShallow                   0.010785004            0.006637748
## Lower Keys\nMesophotic                  0.002930507            0.029957697
## Lower Keys\nShallow                     0.021249108            0.005407182
## Upper Keys\nShallow                     0.019104757            0.014108148
## Upper Keys\nMesophotic                  0.001165877            0.027396216
##                           Riley's Hump\nMesophotic Riley's Hump\nShallow
## Tortugas Bank\nMesophotic                       NA                    NA
## Tortugas Bank\nShallow                          NA                    NA
## Riley's Hump\nMesophotic                        NA                    NA
## Riley's Hump\nShallow                 0.0003676029                    NA
## Lower Keys\nMesophotic                0.0078265244           0.013753336
## Lower Keys\nShallow                   0.0061114733          -0.002212034
## Upper Keys\nShallow                   0.0049598707          -0.001557742
## Upper Keys\nMesophotic                0.0052272451           0.012581289
##                           Lower Keys\nMesophotic Lower Keys\nShallow
## Tortugas Bank\nMesophotic                     NA                  NA
## Tortugas Bank\nShallow                        NA                  NA
## Riley's Hump\nMesophotic                      NA                  NA
## Riley's Hump\nShallow                         NA                  NA
## Lower Keys\nMesophotic                        NA                  NA
## Lower Keys\nShallow                 0.0255855830                  NA
## Upper Keys\nShallow                 0.0211184139        -0.001759316
## Upper Keys\nMesophotic             -0.0002008376         0.024025484
##                           Upper Keys\nShallow Upper Keys\nMesophotic
## Tortugas Bank\nMesophotic                  NA                     NA
## Tortugas Bank\nShallow                     NA                     NA
## Riley's Hump\nMesophotic                   NA                     NA
## Riley's Hump\nShallow                      NA                     NA
## Lower Keys\nMesophotic                     NA                     NA
## Lower Keys\nShallow                        NA                     NA
## Upper Keys\nShallow                        NA                     NA
## Upper Keys\nMesophotic               0.019923                     NA
fk.fst$Pvalues
##                           Tortugas Bank\nMesophotic Tortugas Bank\nShallow
## Tortugas Bank\nMesophotic                        NA                     NA
## Tortugas Bank\nShallow                            0                     NA
## Riley's Hump\nMesophotic                          0                      0
## Riley's Hump\nShallow                             0                      0
## Lower Keys\nMesophotic                            0                      0
## Lower Keys\nShallow                               0                      0
## Upper Keys\nShallow                               0                      0
## Upper Keys\nMesophotic                            0                      0
##                           Riley's Hump\nMesophotic Riley's Hump\nShallow
## Tortugas Bank\nMesophotic                       NA                    NA
## Tortugas Bank\nShallow                          NA                    NA
## Riley's Hump\nMesophotic                        NA                    NA
## Riley's Hump\nShallow                    0.1616162                    NA
## Lower Keys\nMesophotic                   0.0000000                     0
## Lower Keys\nShallow                      0.0000000                     1
## Upper Keys\nShallow                      0.0000000                     1
## Upper Keys\nMesophotic                   0.0000000                     0
##                           Lower Keys\nMesophotic Lower Keys\nShallow
## Tortugas Bank\nMesophotic                     NA                  NA
## Tortugas Bank\nShallow                        NA                  NA
## Riley's Hump\nMesophotic                      NA                  NA
## Riley's Hump\nShallow                         NA                  NA
## Lower Keys\nMesophotic                        NA                  NA
## Lower Keys\nShallow                    0.0000000                  NA
## Upper Keys\nShallow                    0.0000000                   1
## Upper Keys\nMesophotic                 0.8585859                   0
##                           Upper Keys\nShallow Upper Keys\nMesophotic
## Tortugas Bank\nMesophotic                  NA                     NA
## Tortugas Bank\nShallow                     NA                     NA
## Riley's Hump\nMesophotic                   NA                     NA
## Riley's Hump\nShallow                      NA                     NA
## Lower Keys\nMesophotic                     NA                     NA
## Lower Keys\nShallow                        NA                     NA
## Upper Keys\nShallow                        NA                     NA
## Upper Keys\nMesophotic                      0                     NA

Now organize the data to plot

pop.order = c("Riley's Hump\nMesophotic", "Riley's Hump\nShallow", "Tortugas Bank\nMesophotic", "Tortugas Bank\nShallow", "Lower Keys\nMesophotic", "Lower Keys\nShallow", "Upper Keys\nMesophotic", "Upper Keys\nShallow")

# reads in fst matrix
snpFstMa <- as.matrix(fk.fst$Fsts)
upperTriangle(snpFstMa, byrow = TRUE) <- lowerTriangle(snpFstMa)
snpFstMa <- snpFstMa[,pop.order] %>%
  .[pop.order,]
snpFstMa[upper.tri(snpFstMa)] <- NA
snpFstMa <- as.data.frame(snpFstMa)

# rearrange and reformat matrix
snpFstMa$Pop = factor(row.names(snpFstMa), levels = unique(pop.order))

snpQMa <- as.matrix(fk.fst$Pvalues)
upperTriangle(snpQMa, byrow=TRUE) <- lowerTriangle(snpQMa)
snpQMa <- snpQMa[,pop.order] %>%
  .[pop.order,]
snpQMa[upper.tri(snpQMa)] <- NA
snpQMa <- as.data.frame(snpQMa)
snpQMa$Pop = factor(row.names(snpQMa), levels = unique(pop.order))

# melt matrix data (turn from matrix into long dataframe)
snpFstMa$Pop = factor(row.names(snpFstMa), levels = unique(pop.order))
snpFst = melt(snpFstMa, id.vars = "Pop", value.name = "Fst", variable.name = "Pop2", na.rm = FALSE)

snpFst$Fst = round(snpFst$Fst, 3)
snpFst = snpFst %>% mutate(Fst = replace(Fst, Fst < 0, 0))

snpQ = melt(snpQMa, id.vars = "Pop", value.name = "Pval", variable.name = "Pop2", na.rm = FALSE)
# snpQ = snpQ[snpQ$Pop != snpQ$Pop2,]
snpQ$Qval = p.adjust(snpQ$Pval, method = "BH")

snpFst$site = snpFst$Pop
snpFst$site = factor(gsub("\\n.*", "", snpFst$site))
snpFst$site = factor(snpFst$site, levels = levels(snpFst$site)[c(4, 1, 3, 2)])

snpFst$site2 = snpFst$Pop2
snpFst$site2 = factor(gsub("\\n.*", "", snpFst$site2))
snpFst$site2 = factor(snpFst$site2, levels = levels(snpFst$site2)[c(4, 1, 3, 2)])

snpFst$Fst = sprintf('%.3f', snpFst$Fst)
snpFst$Fst = factor(gsub("\\NA", NA, snpFst$Fst))
snpFst$Fst = factor(gsub("\\.000", "", snpFst$Fst))
snpFst$Fst = factor(gsub("\\-", "", snpFst$Fst))

levels(snpFst$Pop) = c("RH\nMeso", "RH\nShallow", "TB\nMeso", "TB\nShallow", "LK\nMeso", "LK\nShallow", "UK\nMeso", "UK \nShallow")
levels(snpFst$Pop2) = c("RH\nMeso", "RH\nShallow", "TB\nMeso", "TB\nShallow", "LK\nMeso", "LK\nShallow", "UK\nMeso", "UK \nShallow")


Plotting fst data as a heatmap

snpHeatmapA = ggplot(data = snpFst, aes(Pop, Pop2, fill = as.numeric(as.character(Fst)))) +
  geom_tile(color = "white") +
  geom_segment(data = snpFst, aes(x = 0.475, xend = -0.5, y = Pop, yend = Pop, color = site), size = 9.75) +
  geom_segment(data = snpFst, aes(x = Pop2, xend = Pop2, y = 0.425, yend = -0.55, color = site2), size = 17) +
  scale_color_manual(values = flPal, guide = NULL) +
  scale_fill_gradient(low = "white", high = "#A14747", limit = c(0, 0.03), space = "Lab", name = expression(paste(italic("F")[ST])), na.value = "white",  guide = "colourbar") +
  geom_text(data = snpFst, aes(Pop, Pop2, label = Fst), color = ifelse(snpQ$Qval <= 0.05,"black", "darkgrey"), size = 3.25, fontface = ifelse(snpQ$Qval < 0.05, "bold", "plain")) +
  guides(fill = guide_colorbar(barwidth = 6, barheight = 0.5, title.position = "top", title.hjust = 0.5)) +
  scale_y_discrete(position = "left") +
  scale_x_discrete(limits = rev(levels(snpFst$Pop))[c(1:8)]) +
  coord_cartesian(xlim = c(1, 8), ylim = c(1, 8), clip = "off") +
  theme_minimal()

snpHeatmap = snpHeatmapA + theme(
  axis.text.x = element_text(vjust = 2, size = 10, hjust = 0.5, color = "gray90"),
  axis.text.y = element_text(size = 10, color = "gray90", margin = margin(r = -5)),
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  panel.grid.major = element_blank(),
  panel.border = element_blank(),
  axis.ticks = element_blank(),
  legend.position = c(0.5, 0.9),
  legend.direction = "horizontal",
  legend.title = element_text(size = 10),
  legend.text = element_text(size = 8)
)

snpHeatmap


ADMIXTURE Plot

Population structure plot from NGSadmix/CLUMPAK results. Based on most probable values for K (2,4).

kColPal = c("#A14747", "#A3764B", "burlywood", "tan3")

admixpops = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone)
admixpops$popdepth = as.factor(paste(admixpops$pop, admixpops$depth, sep = " "))

clumpp2 = read.table("../data/snps/k/ClumppK2.output", header = FALSE)
clumpp4 = read.table("../data/snps/k/ClumppK4.output", header = FALSE)
clumpp2$V1 = admixpops$sample
clumpp4$V1 = admixpops$sample

fkSintAdmix = admixpops %>% left_join(clumpp2[,c(1,6:7)], by = c("sample" = "V1")) %>% 
  left_join(clumpp4[,c(1, 6:9)], by = c("sample" = "V1")) %>% 
  rename("cluster2.1" = "V6.x", "cluster2.2" = "V7.x", "cluster4.1" = "V6.y", "cluster4.2" = "V7.y", "cluster4.3" = "V8", "cluster4.4" = "V9")

fkSintAdmix$pop = factor(fkSintAdmix$pop)
fkSintAdmix$pop = factor(fkSintAdmix$pop, levels(fkSintAdmix$pop)[c( 2, 3, 1, 4)])
fkSintAdmix$depth = factor(fkSintAdmix$depth)
fkSintAdmix$depth = factor(fkSintAdmix$depth, levels(fkSintAdmix$depth)[c(2, 1)])
fkSintAdmix$popdepth = factor(fkSintAdmix$popdepth)
fkSintAdmix$popdepth = factor(fkSintAdmix$popdepth, levels(fkSintAdmix$popdepth)[c(4, 3, 6, 5, 2, 1, 8, 7)])

# fkSintAdmix = arrange(fkSintAdmix, pop, depth, -cluster4.1, cluster4.2, cluster4.3)
fkSintAdmix = arrange(fkSintAdmix, pop, depth, -cluster4.1, -cluster4.2, cluster4.3)
popCounts = fkSintAdmix %>% group_by(popdepth) %>% summarise(n = n())
popCounts
## # A tibble: 8 × 2
##   popdepth                     n
##   <fct>                    <int>
## 1 Riley's Hump Shallow        30
## 2 Riley's Hump Mesophotic     15
## 3 Tortugas Bank Shallow       30
## 4 Tortugas Bank Mesophotic    25
## 5 Lower Keys Shallow          30
## 6 Lower Keys Mesophotic       30
## 7 Upper Keys Shallow          30
## 8 Upper Keys Mesophotic       30
popCountList = reshape2::melt(lapply(popCounts$n,function(x){c(1:x)}))
fkSintAdmix$ord = popCountList$value

fkSintAdmixMelt = melt(fkSintAdmix, id.vars=c("sample", "pop", "depth", "popdepth", "ord"), variable.name="Ancestry", value.name="Fraction")

fkSintAdmixMelt$Ancestry = factor(fkSintAdmixMelt$Ancestry)
fkSintAdmixMelt$Ancestry = factor(fkSintAdmixMelt$Ancestry, levels = rev(levels(fkSintAdmixMelt$Ancestry)))

popAnno = data.frame(x1 = c(0.25, 0.25, 0.25, 0.25), x2 = c(30.75, 30.75, 30.75, 30.75),
                     y1 = -0.14, y2 = -0.14, sample = NA, Ancestry = NA, depth = "Mesophotic", 
                     ord  = NA, Fraction = NA,
                     pop = c("Riley's Hump", "Tortugas Bank", 
                                  "Lower Keys", "Upper Keys"))
popAnno$pop = factor(popAnno$pop)
popAnno$pop = factor(popAnno$pop, levels = levels(popAnno$pop)[c(4, 1, 3, 2)])

admixTheme = theme_bw()+
  theme(plot.title = element_text(hjust = 0.5, size = 10),
  panel.grid = element_blank(),
  panel.background = element_rect(fill = "gray25", colour = "gray25"),
  panel.border = element_rect(fill = NA, color = "black", size = 1, linetype = "solid"),
  panel.spacing.x = grid:::unit(0.05, "lines"),
  panel.spacing.y = grid:::unit(0.05, "lines"),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_blank(),
  axis.ticks.y = element_blank(),
  axis.title = element_blank(),
  strip.background.x = element_blank(),
  strip.background.y = element_blank(),
  strip.text = element_text(size = 10),
  strip.text.y.left = element_text(size = 10, angle = 90),
  strip.text.x.bottom = element_text(vjust = -.1, color = "grey90"),
  legend.key = element_blank(),
  legend.position = "none",
  legend.title = element_text(size = 8))

admixPlotK2A = ggplot(data = subset(fkSintAdmixMelt, subset = fkSintAdmixMelt$Ancestry %in% c("cluster2.1", "cluster2.2")), aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
    geom_segment(data = popAnno, aes(x = x1, xend = x2, y = y1, yend = y2, color = pop), size = 8) +
  geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
  facet_grid(factor(depth) ~ pop, scales = "free", switch = "both", space = "free") +
  scale_x_discrete(expand = c(0.01, 0.01)) +
  scale_y_continuous(expand = c(0.005, 0.005)) +
  scale_fill_manual(values = kColPal[c(1,3)]) +
  scale_color_manual(values = flPal) +
  # ggtitle(expression(paste(italic("K")," = 2", sep =""))) +
  # coord_cartesian(expand = TRUE, clip = "off")
  coord_cartesian(ylim = c(-0.01, 1.01), clip = "off")

admixPlotK2 = admixPlotK2A + admixTheme

admixPlotK4A = ggplot(data = subset(fkSintAdmixMelt, subset = !(fkSintAdmixMelt$Ancestry %in% c("cluster2.1", "cluster2.2"))), aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
    geom_segment(data = popAnno, aes(x = x1, xend = x2, y = y1, yend = y2, color = pop), size = 8) +
  geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
  facet_grid(factor(depth) ~ pop, scales = "free", switch = "both", space = "free") +
  scale_x_discrete(expand = c(0.01, 0.01)) +
  scale_y_continuous(expand = c(0.005, 0.005)) +
  scale_fill_manual(values = kColPal) +
  scale_color_manual(values = flPal) +
  # ggtitle(expression(paste(italic("K")," = 4", sep =""))) +
  # coord_cartesian(expand = TRUE, clip = "off")
  coord_cartesian(ylim = c(-.01, 1.01), clip = "off")
 
admixPlotK4 = admixPlotK4A + admixTheme
# admixPlotK4 = admixPlotK4A + admixTheme + 
#   theme(strip.text.y.left = element_blank())

admixPlot = (admixPlotK2 / admixPlotK4)

admixPlot

# pie

# fkSintAdmixMelt %>% filter(depth == "Shallow", Ancestry != "cluster2.1", Ancestry != "cluster2.2") %>% group_by(depth, Ancestry) %>% summarise(Fraction = sum(Fraction)/n()) %>%
# ggplot(data = ., aes(x = depth, y = Fraction, fill = Ancestry)) +
#   geom_bar(stat="identity", width=1) +
#   coord_polar("y", start=0) +
#   scale_fill_manual(values = rev(kColPal)) +
#   theme_void()

# fkSintAdmixMelt %>% filter(depth == "Mesophotic", Ancestry != "cluster2.1", Ancestry != "cluster2.2") %>% group_by(depth, Ancestry) %>% summarise(Fraction = sum(Fraction)/n()) %>%
# ggplot(data = ., aes(x = depth, y = Fraction, fill = Ancestry)) +
#   geom_bar(stat="identity", width=1) +
#   coord_polar("y", start=0) +
#   scale_fill_manual(values = rev(kColPal)) +
#   theme_void()

Population structure plots

sintStructureFig = ((((sintDapcPlot)+ theme(axis.title.y = element_text(margin = margin(r = -45, unit = "pt")), legend.text = element_text(size = 10), legend.title = element_text(size = 10), legend.spacing = unit(2, "pt"),legend.key.size = unit(15, "pt"), legend.background = element_blank())) | ((sintDapcDepthPlot) + theme(legend.text = element_text(size = 10), legend.title = element_text(size = 10), legend.spacing = unit(2, "pt"),legend.key.size = unit(15, "pt"), legend.background = element_blank()))) / ((snpHeatmap + theme(axis.text.x = element_text(vjust = 11.5, size = 10, hjust = 0.5, color = "gray90"), axis.text.y = element_text(size = 10, color = "gray90", margin = margin(r = -5, unit = "pt")))) | (admixPlotK4))) + plot_layout(heights = c(1, 0.8)) + plot_annotation(tag_levels = "A")

# sintStructureFig

ggsave("../figures/figure4.png", plot = sintStructureFig, height = 18, width = 28, units = "cm", dpi = 300)
ggsave("../figures/figure4.svg", plot = sintStructureFig, height = 18, width = 28, units = "cm", dpi = 300)

Isolation by distance and redundancy analysis

S. intersepta isolation by distance

Procrustes analysis using IBS distance matrix and geographic locations

sintIBS = read.table("../data/snps/sintNoClones.ibsMat")%>% as.matrix() %>% as.dist(diag = FALSE)

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "Site" = site, "Depth" = depthZone) %>% mutate(popDepth = paste(Site, Depth))

coords = fknmsSites %>% group_by(site, depthZone, siteID) %>% summarise(latDD = first(latDD), longDD = first(longDD)) %>% filter(siteID %in% c("Ian's Lumps Site 52", "Site 48", "Site 47", "Site 45", "Site 35/36", "Site 37", "Site 39", "Site 19")) %>% dplyr::select(-site) %>% mutate(popDepth = paste(site, depthZone)) %>% droplevels()
## `summarise()` has grouped output by 'site', 'depthZone'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `site`
coords = coords[-5,]

sintIBD = popData %>% left_join(coords) %>% dplyr::select(-site, -depthZone, -popDepth)
## Joining, by = "popDepth"
rownames(sintIBD) = sintIBD$sample

sintProj = st_as_sf(x = sintIBD, coords = c("longDD", "latDD"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>% st_transform(crs = 5070) 

sintIBD2 = sintIBD %>% dplyr::select(-longDD, -latDD) %>% cbind(st_coordinates(sintProj))

set.seed(694) 
ibd = protest(X = sintIBD2[,c(5, 6)], Y = cmdscale(sintIBS), choices = c(1, 2), permutations = 9999, symmetric = FALSE)

ibd
## 
## Call:
## protest(X = sintIBD2[, c(5, 6)], Y = cmdscale(sintIBS), permutations = 9999,      choices = c(1, 2), symmetric = FALSE) 
## 
## Procrustes Sum of Squares (m12 squared):        0.9893 
## Correlation in a symmetric Procrustes rotation: 0.1034 
## Significance:  0.1267 
## 
## Permutation: free
## Number of permutations: 9999
plot(ibd, kind = 2)

plot(ibd, kind = 1)

We can make a plot demonstrating the lack of IBD from the procruestes analysis

ibdPlot = procrustes(X = sintIBD2[,c(5, 6)], Y = cmdscale(sintIBS), choices = c(1, 2), symmetric = FALSE)

sintIBDPlotData = cbind(ibdPlot$X, ibdPlot$Yrot) %>% as.data.frame()
rownames(sintIBDPlotData) = rownames(ibdPlot$X)
colnames(sintIBDPlotData) = c("lon", "lat", "ibs1", "ibs2")
sintIBDPlotData$sample = row.names(sintIBDPlotData)
sintIBDPlotData$sample = gsub(pattern = "\\.2", "", sintIBDPlotData$sample)
sintIBDPlotData = sintIBDPlotData %>% left_join(popData) %>% relocate(sample, .before = lon)
## Joining, by = "sample"
sintIBDPlotData$Depth = factor(sintIBDPlotData$Depth)
sintIBDPlotData$Depth = factor(sintIBDPlotData$Depth, levels(sintIBDPlotData$Depth)[c(2,1)])
sintIBDPlotData$Site = factor(sintIBDPlotData$Site)
sintIBDPlotData$Site = factor(sintIBDPlotData$Site, levels(sintIBDPlotData$Site)[c(4, 1, 3, 2)])
  
head(sintIBDPlotData)
##   sample       lon       lat       ibs1     ibs2          Site      Depth
## 1 SFK001 -120448.3 -27548.28   5841.166 4337.123 Tortugas Bank Mesophotic
## 2 SFK002 -120448.3 -27548.28 -13609.418 7171.627 Tortugas Bank Mesophotic
## 3 SFK003 -120448.3 -27548.28 -13679.616 7422.906 Tortugas Bank Mesophotic
## 4 SFK004 -120448.3 -27548.28   5020.850 5091.913 Tortugas Bank Mesophotic
## 5 SFK005 -120448.3 -27548.28   5276.879 5081.843 Tortugas Bank Mesophotic
## 6 SFK006 -120448.3 -27548.28   4718.687 5710.928 Tortugas Bank Mesophotic
##                   popDepth
## 1 Tortugas Bank Mesophotic
## 2 Tortugas Bank Mesophotic
## 3 Tortugas Bank Mesophotic
## 4 Tortugas Bank Mesophotic
## 5 Tortugas Bank Mesophotic
## 6 Tortugas Bank Mesophotic
convertedIBD = sintIBDPlotData %>% mutate(lon = lon+ibdPlot$translation[1], lat = lat+ibdPlot$translation[2], ibs1 = ibs1+ibdPlot$translation[1], ibs2 = ibs2+ibdPlot$translation[2])

sintIBDPlotPops = convertedIBD %>% dplyr::select(Site, Depth, lon, lat) %>% group_by(lon, lat)

#Calculate the angle of rotation, and then the slope of the rotated axis
ibdTheta = acos(ibdPlot$rotation[1,1]) 
ibdSlope = tan(ibdTheta)

#Calculate the y-intercepts for rotated axes
ibdInt = ibdPlot$translation[2] - (ibdSlope*ibdPlot$translation[1])
ibdInt2 = ibdPlot$translation[2] - (-1/ibdSlope*ibdPlot$translation[1])

florida2 = st_transform(florida, crs = 5070)

sintIBDPlotA = ggplot() +
  geom_sf(data = florida2, fill = "white", size = 0.25) +
  geom_abline(intercept = ibdInt, slope = ibdSlope, color = "gray75", linetype = 2) +
  geom_abline(intercept = ibdInt2, slope = -(1/ibdSlope), color = "gray75", linetype = 2) +
    geom_segment(data = convertedIBD, aes(x = ibs1, y = ibs2, xend = lon, yend = lat, color = Site), alpha = 0.5) +
  geom_point(data = convertedIBD, aes(x = ibs1, y = ibs2, color = Site), shape = 16, alpha = 0.5)+
  geom_point(data = sintIBDPlotPops, aes(x = lon, y = lat, fill = Site, shape = Depth), size = 2) +
    annotate(geom = "label", x = 1344053.08266 , y = 368848.804494, label = "           ", size = 10) +
  annotate(geom = "text", x = 1340417.86738 , y = 372764.041467, label = "Procrustes analysis:", size = 3) +
  annotate(geom = "text", x = 1345669.60498, y = 364652.694261, label = "italic(t[0]) == 0.0678*','~italic(p) == 0.5173", parse = TRUE, size = 3) +
  scale_color_manual(values = flPal) +
  scale_fill_manual(values = flPal) +
  scale_shape_manual(values = c(24, 25), name = "Depth zone") +
  coord_sf(default_crs = sf::st_crs(5070), ylim = c(250000, 375000), xlim = c(1320000, 1600000)) +
  annotation_scale(location = "br") +
  annotation_north_arrow(location = "br", which_north = "true", style = north_arrow_minimal(), pad_x = unit(-0.25, "cm") , pad_y = unit(0.75, "cm")) +
  guides(color = "none", fill = guide_legend(override.aes = list(shape = 15, color = flPal, size = 3), ncol =2, order = 1), shape = guide_legend(ncol = 1)) +
  theme_bw()

sintIBDPlot = sintIBDPlotA +
  theme(panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        axis.title = element_blank(),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(color = "black", size = 8),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.box = "horizontal",
        legend.key = element_blank(),
        legend.background = element_blank(),
        legend.title = element_text(color = "black", size = 8),
        legend.text = element_text(color = "black", size = 8)
        )

sintIBDPlot

ggsave("../figures/figureS4.png", plot = sintIBDPlot, height = 10, width = 16, units = "cm", dpi = 300)
ggsave("../figures/figureS4.svg", plot = sintIBDPlot, height = 8, width = 14, units = "cm", dpi = 300)


Distance-based redundancy analysis (dbRDA)

# Set directory for sdmPredictors to download/search for previously downloaded datasets 
# Also allow more than 60 seconds to download larger datasets before timing out
options(sdmpredictors_datadir="../data/snps/bioOracle", timeout = max(300, getOption("timeout")))

sintMa = as.dist(read.table("../data/snps/sintNoClones.ibsMat"))
popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select(Sample, site, depthZone, latDD, longDD, depthM)

datasets = list_datasets(terrestrial = FALSE, marine = TRUE, freshwater = FALSE)

# Extract present-day data sets
present = list_layers(datasets) %>% dplyr::select(dataset_code, layer_code, name, units, description, contains("cellsize"), version) %>% filter(version == 2.2) %>% filter(layer_code %in% c("BO22_damean", "BO22_parmean", "BO22_ph", "BO22_curvelmax_bdmean", "BO22_salinitymean_bdmean", "BO22_salinitymean_ss", "BO22_curvelmean_ss", "BO22_curvelmean_bdmean", "BO22_dissoxmean_bdmean", "BO22_lightbotmax_bdmean", "BO22_lightbotmean_bdmean", "BO22_nitratemean_bdmean", "BO22_tempmax_bdmean", "BO22_tempmean_bdmean", "BO22_tempmean_ss", "BO22_ppmean_ss", "BO22_dissoxmean_ss", "BO22_ppmean_bdmean", "BO22_chlomean_ss", "BO22_chlomean_bdmean"))

envVar = load_layers(present$layer_code)

envData = data.frame(popData, raster::extract(envVar, popData[,5:4]))

# See what data are highly colinear
corData = rcorr(as.matrix(envData[,c(6:ncol(envData))]), type = "pearson")
corDataFlat = melt(corData$r, value.name = "r")
pDataFlat = melt(corData$P, value.name = "p")
corDataBind = corDataFlat %>% left_join(pDataFlat, by = c("Var1","Var2"))

ggplot(corDataBind) +
  geom_tile(aes(x = Var1, y = Var2, fill = r)) +
  scale_fill_gradient2(low = "#3B9AB2FF", high = "#F21A00FF") +
  geom_text(data = filter(corDataBind, r >= 0.7, p < 0.05),aes(x = Var1, y = Var2, label = round(r, 2))) +
    geom_text(data = filter(corDataBind, r <= -0.7, p < 0.05),aes(x = Var1, y = Var2, label = round(r, 2))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

# Removing correlated data and biologically irrelevant measurements
envData2 = envData %>% dplyr::select("BO22_curvelmean_ss", "depthM", "BO22_lightbotmean_bdmean", "BO22_tempmean_bdmean", "BO22_ppmean_bdmean", "BO22_dissoxmean_bdmean")

# Checking again to make sure we've eliminated all colinearity
corData2 = cor(envData2)
corMelt2 = melt(corData2)

ggplot(corMelt2) +
  geom_tile(aes(x = Var1, y = Var2, fill = value)) +
  scale_fill_gradient2(low = "#3B9AB2FF", high = "#F21A00FF") +
  geom_text(data = corMelt2[corMelt2$value >= 0.7,],aes(x = Var1, y = Var2, label = round(value, 2))) +
    geom_text(data = corMelt2[corMelt2$value <= -0.7,],aes(x = Var1, y = Var2, label = round(value, 2))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

Looks good now!

Now we check variance inflation factors(VIF) to see if there is any multi-colinearity (VIF ≥ 10) between explanatory variables

vif = diag(solve(cor(envData2)))
vif
##       BO22_curvelmean_ss                   depthM BO22_lightbotmean_bdmean 
##                 4.040956                 1.694737                 8.342723 
##     BO22_tempmean_bdmean       BO22_ppmean_bdmean   BO22_dissoxmean_bdmean 
##                 6.946536                 8.158057                 6.204014

No issues with the VIF measures so we can continue to run the redundancy analysis

dbrda0 = dbrda(sintMa ~ 1, data = envData2)
dbrda1 = dbrda(sintMa ~ ., data = envData2)
anova(dbrda1)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = sintMa ~ BO22_curvelmean_ss + depthM + BO22_lightbotmean_bdmean + BO22_tempmean_bdmean + BO22_ppmean_bdmean + BO22_dissoxmean_bdmean, data = envData2)
##           Df SumOfSqs      F Pr(>F)    
## Model      6   0.3286 1.8175  0.001 ***
## Residual 213   6.4178                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(dbrda1)
## $r.squared
## [1] 0.04870246
## 
## $adj.r.squared
## [1] 0.02190535

The model explains some of the genetic variation we see. Now we can run forward selection to find the best model that explains the data without over-fitting.

set.seed(092)
bestDbrda = ordiR2step(dbrda0, dbrda1, permutations = how(nperm = 9999), R2permutations = 9999, direction = "forward")
## Step: R2.adj= 0 
## Call: sintMa ~ 1 
##  
##                                R2.adjusted
## <All variables>             0.021905350693
## + depthM                    0.016830338892
## + BO22_lightbotmean_bdmean  0.002294998071
## + BO22_dissoxmean_bdmean    0.001700807951
## + BO22_ppmean_bdmean        0.001687667218
## + BO22_tempmean_bdmean      0.000004163804
## <none>                      0.000000000000
## + BO22_curvelmean_ss       -0.000254882481
## 
##          Df    AIC      F Pr(>F)    
## + depthM  1 418.24 4.7489 0.0001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.01683034 
## Call: sintMa ~ depthM 
##  
##                            R2.adjusted
## <All variables>             0.02190535
## + BO22_dissoxmean_bdmean    0.01958483
## + BO22_ppmean_bdmean        0.01897755
## + BO22_tempmean_bdmean      0.01839357
## + BO22_curvelmean_ss        0.01687645
## + BO22_lightbotmean_bdmean  0.01685398
## <none>                      0.01683034
## 
##                          Df    AIC      F Pr(>F)   
## + BO22_dissoxmean_bdmean  1 418.61 1.6125 0.0023 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.01958483 
## Call: sintMa ~ depthM + BO22_dissoxmean_bdmean 
##  
##                            R2.adjusted
## <All variables>             0.02190535
## + BO22_tempmean_bdmean      0.01996180
## + BO22_ppmean_bdmean        0.01985213
## + BO22_curvelmean_ss        0.01966802
## <none>                      0.01958483
## + BO22_lightbotmean_bdmean  0.01899736
## 
##                        Df    AIC      F Pr(>F)
## + BO22_tempmean_bdmean  1 419.51 1.0835 0.2004
bestDbrda$anova
##                            R2.adj Df    AIC      F Pr(>F)    
## + depthM                 0.016830  1 418.24 4.7489 0.0001 ***
## + BO22_dissoxmean_bdmean 0.019585  1 418.61 1.6125 0.0023 ** 
## <All variables>          0.021905                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Depth and Primary productivity are the environmental variables that best explain the variation (although only 2.18%).

model = envData2 %>% dplyr::select("Depth" = depthM, "DO" = BO22_dissoxmean_bdmean)

sintDbrda = dbrda(sintMa ~ Depth + DO, model)

dbrdaVarPart = varpart(as.dist(sintMa), ~ Depth, ~ DO, data = model)
dbrdaDepth = dbrda(sintMa ~ Depth, data = model)
dbrdaDO = dbrda(sintMa ~ DO, data = model)

dbrdaVarPart
## 
## Partition of squared Unknown user-supplied distance in dbRDA 
## 
## Call: varpart(Y = as.dist(sintMa), X = ~Depth, ~DO, data = model)
## 
## Explanatory tables:
## X1:  ~Depth
## X2:  ~DO 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 6.7463 
## No. of observations: 220 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+b] = X1            1   0.02132       0.01683     TRUE
## [b+c] = X2            1   0.00626       0.00170     TRUE
## [a+b+c] = X1+X2       2   0.02854       0.01958     TRUE
## Individual fractions                                    
## [a] = X1|X2           1                 0.01788     TRUE
## [b]                   0                -0.00105    FALSE
## [c] = X2|X1           1                 0.00275     TRUE
## [d] = Residuals                         0.98042    FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
set.seed(003)
anova(sintDbrda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = sintMa ~ Depth + DO, data = model)
##           Df SumOfSqs      F Pr(>F)    
## Model      2   0.1925 3.1874  0.001 ***
## Residual 217   6.5538                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(002)
anova.cca(dbrdaDepth)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = sintMa ~ Depth, data = model)
##           Df SumOfSqs      F Pr(>F)    
## Model      1   0.1438 4.7489  0.001 ***
## Residual 218   6.6025                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(001)
anova(dbrdaDO)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = sintMa ~ DO, data = model)
##           Df SumOfSqs      F Pr(>F)  
## Model      1   0.0422 1.3731  0.025 *
## Residual 218   6.7041                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now create the best model and prepare to plot

sintRdaVar = round(sintDbrda$CA$eig/sum(sintDbrda$CA$eig)*100, 1)
head(sintRdaVar)
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 
##  8.1  3.9  3.3  0.7  0.7  0.7
sintRdaVarFit = round(sintDbrda$CCA$eig/sum(sintDbrda$CCA$eig)*100, 1)
head(sintRdaVarFit)
## dbRDA1 dbRDA2 
##     79     21
sintI2P = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66,68,164,166,209,211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone)

sintI2P$popdepth = paste(sintI2P$pop, sintI2P$depth)

sintRdaPoints = as.data.frame(scores(sintDbrda)$sites)
sintRdaPoints$sample = sintI2P$sample
head(sintRdaPoints)
##        dbRDA1     dbRDA2 sample
## V1 -0.5589827 -0.3199523 SFK001
## V2 -0.1680224 -0.7858262 SFK002
## V3 -0.1728142 -0.8276181 SFK003
## V4 -0.5436297 -0.3341232 SFK004
## V5 -0.5799407 -0.2517953 SFK005
## V6 -0.5882035 -0.2476965 SFK006
sintDbrdaData1 = sintI2P %>% left_join(sintRdaPoints) %>% cbind(K = dapc1$grp)
## Joining, by = "sample"
# sintDbrdaData1 = sintI2P %>% left_join(sintRdaPoints)
head(sintDbrdaData1)
##                         sample           pop      depth                 popdepth
## fk_S001.trim.un.bt2.bam SFK001 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
## fk_S002.trim.un.bt2.bam SFK002 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
## fk_S003.trim.un.bt2.bam SFK003 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
## fk_S004.trim.un.bt2.bam SFK004 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
## fk_S005.trim.un.bt2.bam SFK005 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
## fk_S006.trim.un.bt2.bam SFK006 Tortugas Bank Mesophotic Tortugas Bank Mesophotic
##                             dbRDA1     dbRDA2 K
## fk_S001.trim.un.bt2.bam -0.5589827 -0.3199523 1
## fk_S002.trim.un.bt2.bam -0.1680224 -0.7858262 2
## fk_S003.trim.un.bt2.bam -0.1728142 -0.8276181 2
## fk_S004.trim.un.bt2.bam -0.5436297 -0.3341232 1
## fk_S005.trim.un.bt2.bam -0.5799407 -0.2517953 1
## fk_S006.trim.un.bt2.bam -0.5882035 -0.2476965 1
tail(sintDbrdaData1)
##                         sample        pop      depth              popdepth     dbRDA1
## fk_S215.trim.un.bt2.bam SFK215 Upper Keys Mesophotic Upper Keys Mesophotic -0.5551502
## fk_S216.trim.un.bt2.bam SFK216 Upper Keys Mesophotic Upper Keys Mesophotic -0.3906943
## fk_S217.trim.un.bt2.bam SFK217 Upper Keys    Shallow    Upper Keys Shallow -0.4665246
## fk_S218.trim.un.bt2.bam SFK218 Upper Keys    Shallow    Upper Keys Shallow  1.5851168
## fk_S219.trim.un.bt2.bam SFK219 Upper Keys    Shallow    Upper Keys Shallow -0.4098178
## fk_S220.trim.un.bt2.bam SFK220 Upper Keys    Shallow    Upper Keys Shallow -0.4142953
##                            dbRDA2 K
## fk_S215.trim.un.bt2.bam 0.4903874 1
## fk_S216.trim.un.bt2.bam 0.2072862 1
## fk_S217.trim.un.bt2.bam 0.6889664 1
## fk_S218.trim.un.bt2.bam 1.3642264 4
## fk_S219.trim.un.bt2.bam 0.5406104 1
## fk_S220.trim.un.bt2.bam 0.7035300 1
envLoad = as.data.frame(sintDbrda$CCA$biplot)
envLoad$var = row.names(envLoad)

sintDbrdaData = merge(sintDbrdaData1, aggregate(cbind(mean.x = dbRDA1, mean.y = dbRDA2)~popdepth, sintDbrdaData1, mean), by="popdepth") %>% merge(., aggregate(cbind(mean.x.K = dbRDA1, mean.y.K = dbRDA2)~K, sintDbrdaData1, mean), by="K")
# sintDbrdaData = merge(sintDbrdaData1, aggregate(cbind(mean.x = dbRDA1, mean.y = dbRDA2)~popdepth, sintDbrdaData1, mean), by="popdepth")

sintDbrdaData$depth = factor(sintDbrdaData$depth)
sintDbrdaData$depth = factor(sintDbrdaData$depth, levels(sintDbrdaData$depth)[c(2,1)])
sintDbrdaData$pop = factor(sintDbrdaData$pop)
sintDbrdaData$pop = factor(sintDbrdaData$pop, levels(sintDbrdaData$pop)[c(4, 1, 3, 2)])

head(sintDbrdaData)
##   K              popdepth sample        pop      depth     dbRDA1      dbRDA2
## 1 1 Lower Keys Mesophotic SFK101 Lower Keys Mesophotic -0.5922790 -0.14670222
## 2 1 Lower Keys Mesophotic SFK102 Lower Keys Mesophotic -0.4702825 -0.13777606
## 3 1 Lower Keys Mesophotic SFK103 Lower Keys Mesophotic -0.4623982 -0.14996794
## 4 1 Lower Keys Mesophotic SFK104 Lower Keys Mesophotic -0.5617873 -0.02281652
## 5 1 Lower Keys Mesophotic SFK105 Lower Keys Mesophotic -0.5706488 -0.05593455
## 6 1 Lower Keys Mesophotic SFK106 Lower Keys Mesophotic -0.5169524 -0.11334307
##       mean.x     mean.y   mean.x.K   mean.y.K
## 1 -0.5036297 -0.1133085 -0.4813059 0.08451714
## 2 -0.5036297 -0.1133085 -0.4813059 0.08451714
## 3 -0.5036297 -0.1133085 -0.4813059 0.08451714
## 4 -0.5036297 -0.1133085 -0.4813059 0.08451714
## 5 -0.5036297 -0.1133085 -0.4813059 0.08451714
## 6 -0.5036297 -0.1133085 -0.4813059 0.08451714
domClust = fkSintAdmix %>% dplyr::select(-"cluster2.1", -"cluster2.2", -"ord", -"pop", -"depth", -"popdepth") %>% left_join(het) %>% mutate(cluster = ifelse(cluster4.1 < 0.75 & cluster4.2 < 0.75  & cluster4.3 < 0.75 & cluster4.4 < 0.75, "NA", ifelse(cluster4.1 >=0.75, 1, ifelse(cluster4.2 >= 0.75, 2, ifelse(cluster4.3 >= 0.75, 3,ifelse(cluster4.4 >= 0.75, 4, 0))))))
## Joining, by = "sample"
domClust$cluster = as.factor(domClust$cluster)
levels(domClust$cluster) = c("Red", "Brown", "Tan", "Orange", "Admixed")

sintDbrdaData2 = domClust %>% dplyr::select(sample, cluster) %>% left_join(sintDbrdaData)
## Joining, by = "sample"

Now plot the dbRDA

flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]
kColPal = c("#A14747", "#A3764B", "burlywood", "tan3")

sintDbrdaPlotA = ggplot() +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  geom_point(data = sintDbrdaData2, aes(x = dbRDA1, y = dbRDA2, fill = pop, shape = depth, color = pop), size = 2, alpha = 0.5, show.legend = FALSE) +
  scale_shape_manual(values = c(24,25), name = "Depth Zone") +
  geom_segment(data = envLoad, aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2), color = "gray35", arrow = arrow(length = unit(0.15, "cm"), type = "open"), size = 0.65) +
  geom_text(data = envLoad[1,], aes(x = dbRDA1+0.05, y = dbRDA2+0.2, label = var), color = "gray35", size = 3) +
   geom_text(data = envLoad[2,], aes(x = dbRDA1+0.1, y = dbRDA2, label = var), color = "gray35", size = 3) +
  geom_point(data = sintDbrdaData2, aes(x = mean.x, y = mean.y, fill = pop, shape = depth), size = 3, color = "black") + #population centroids indicated by triangles
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = flPal, name = "Site") +
  labs(x = paste0("dbRDA 1 (", sintRdaVarFit[1], "% [",sintRdaVar[1], "%])"), y = paste0("dbRDA 2 (", sintRdaVarFit[2], "% [",sintRdaVar[2], "%])")) +
  guides(shape = guide_legend(override.aes = list(size = 2, stroke = 0.5, alpha = NA), order = 2), fill = guide_legend(override.aes = list(shape = 22, size = 2, color = flPal, alpha = NA), order = 1))+
  theme_bw()

sintDbrdaPlot1 = sintDbrdaPlotA +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.position = "right",
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

sintDbrdaPlotB = ggplot() +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  geom_point(data = sintDbrdaData2, aes(x = dbRDA1, y = dbRDA2, fill = cluster, shape = depth, color = cluster), size = 2, alpha = 0.5) +
  scale_shape_manual(values = c(24, 25), name = "Depth Zone") +
  geom_segment(data = envLoad, aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2), color = "gray35", arrow = arrow(length = unit(0.15, "cm"), type = "open"), size = 0.65) +
  geom_text(data = envLoad[1,], aes(x = dbRDA1+0.05, y = dbRDA2+0.2, label = var), color = "gray35", size = 3) +
   geom_text(data = envLoad[2,], aes(x = dbRDA1+0.1, y = dbRDA2, label = var), color = "gray35", size = 3) +
  scale_fill_manual(values = c(kColPal, "gray50"), name = "Lineage") +
  scale_color_manual(values = c(kColPal, "gray50"), name = "Lineage", guide = NULL) +
  labs(x = paste0("dbRDA 1 (", sintRdaVarFit[1], "% [",sintRdaVar[1], "%])"), y = paste0("dbRDA 2 (", sintRdaVarFit[2], "% [",sintRdaVar[2], "%])")) +
  guides(shape = guide_legend(override.aes = list(size = 2, stroke = 0.5, alpha = 1), order = 2), fill = guide_legend(override.aes = list(shape = 22, size = 2, alpha = NA, color = c(kColPal, "gray50")), order = 1))+
  theme_bw()

sintDbrdaPlot2 = sintDbrdaPlotB +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.position = "right",
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

dbrdaPlots = (sintDbrdaPlot1 / sintDbrdaPlot2) +
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A')& 
  theme(legend.spacing = unit(2, "pt"),
        legend.key.height = unit(10, "pt"),
        plot.tag = element_text(color = "black"))

dbrdaPlots

ggsave("../figures/figure5.png", plot = dbrdaPlots, height = 14, width = 14, units = "cm", dpi = 300)
ggsave("../figures/figure5.svg", plot = dbrdaPlots, height = 14, width = 14, units = "cm", dpi = 300)


Genetic connectivity

Checking deviance among model runs from BayesAss we ran on HPC

fileList = substr(list.files("../data/snps/bayesAss/", "BA3trace.*.txt$"),1,10)

bayesian_deviance <- function(trace, burnin = 0, sampling.interval = 0){
  if(burnin == 0) stop('No burnin specified')
  if(sampling.interval == 0) stop('No sampling interval specified')
  range <- (trace$State > burnin & trace$State %% sampling.interval == 0)
  D <- -2*mean(trace$LogProb[range])
  return(D)
}

for(i in 1:length(fileList)){
  assign(fileList[i], read.delim(paste("../data/snps/bayesAss/", fileList[i], ".txt", sep = ""))) %>% dplyr::select(-last_col())
  print(paste(fileList[i], bayesian_deviance(get(fileList[i]), burnin = 3000000, sampling.interval = 100)))
}
## [1] "BA3trace01 6276550.568"
## [1] "BA3trace02 6274555.94971429"
## [1] "BA3trace03 6273601.16057143"
## [1] "BA3trace04 6274629.62228571"
## [1] "BA3trace05 6273869.57542857"
## [1] "BA3trace06 6273118.99771429"
## [1] "BA3trace07 6274232.18228571"
## [1] "BA3trace08 6274603.30514286"
## [1] "BA3trace09 6273023.81885714"
## [1] "BA3trace10 6273319.95257143"

All traces have similar deviance (this is good!). Using the trace with the lowest deviance (BA3trace09.txt, in this case)

bayesAss = read.delim("../data/snps/bayesAss/BA3trace09.txt") %>% filter(State > 3000000) %>% dplyr::select(-State, -LogProb, -X)

baMean = bayesAss %>% summarise(across(everything(), list(mean))) %>% t() %>% as_tibble() %>% rename(., mean=V1) %>% mutate(pops = colnames(bayesAss))

baSumm = bayesAss %>% summarise(across(everything(), list(median))) %>% t() %>% as.data.frame() %>% rename(., median=V1) %>% mutate(pops = baMean$pops, mean = round(baMean$mean, 3)) %>% relocate(median, .after = mean)

baSumm$median = round(baSumm$median, 3)

baHpd =as.data.frame(t(sapply(bayesAss, emp.hpd)))
colnames(baHpd) = c("hpdLow", "hpdHigh")
baHpd$pops = rownames(baHpd)

ESS = as.data.frame(sapply(bayesAss, ESS))
colnames(ESS) = "ESS"

baSumm = baSumm %>% left_join(baHpd)
## Joining, by = "pops"
baSumm$hpdLow = round(baSumm$hpdLow, 3)
baSumm$hpdHigh = round(baSumm$hpdHigh, 3)
baSumm$ESS = ESS$ESS

### FROM BAYESASS: ###
## Population Index -> Population Label:
## 0->TortugasBank_Mesophotic 1->TortugasBank_Shallow
## 2->RileysHump_Mesophotic 3->RileysHump_Shallow
## 4->LowerKeys_Mesophotic 5->LowerKeys_Shallow
## 6->UpperKeys_Shallow 7->UpperKeys_Mesophotic

popi = rep(c("Tortugas Bank\nMesophotic", "Tortugas Bank\nShallow", "Riley's Hump\nMesophotic", "Riley's Hump\nShallow", "Lower Keys\nMesophotic", "Lower Keys\nShallow", "Upper Keys\nShallow", "Upper Keys\nMesophotic"), each = 8)

popj = rep(c("Tortugas Bank\nMesophotic", "Tortugas Bank\nShallow", "Riley's Hump\nMesophotic", "Riley's Hump\nShallow", "Lower Keys\nMesophotic", "Lower Keys\nShallow", "Upper Keys\nShallow", "Upper Keys\nMesophotic"), times = 8)

baSumm = baSumm %>% mutate(pop.i = popi, pop.j = popj) %>% relocate(c(pop.i, pop.j), .after = pops) %>% dplyr::select(-pops)

baSumm$pop.i = factor(baSumm$pop.i)
baSumm$pop.i = factor(baSumm$pop.i, levels = levels(baSumm$pop.i)[c(8, 2, 6, 4, 7, 1, 5, 3)])

baSumm$pop.j = factor(baSumm$pop.j)
baSumm$pop.j = factor(baSumm$pop.j, levels = levels(baSumm$pop.j)[c(8, 2, 6, 4, 7, 1, 5, 3)])

baSumm$site.i = word(baSumm$pop.i, 1, sep = "\n")
baSumm$site.i = factor(baSumm$site.i)
baSumm$site.i = factor(baSumm$site.i, levels = levels(baSumm$site.i)[c(4, 1, 3, 2)])

baSumm$site.j = word(baSumm$pop.j, 1, sep = "\n")
baSumm$site.j = factor(baSumm$site.j)
baSumm$site.j = factor(baSumm$site.j, levels = levels(baSumm$site.j)[c(4, 1, 3, 2)])

baSumm$depth.i = word(baSumm$pop.i, 2, sep = "\n")
baSumm$depth.i = factor(baSumm$depth.i)
baSumm$depth.i = factor(baSumm$depth.i, levels = levels(baSumm$depth.i)[c(2, 1)])

baSumm$depth.j = word(baSumm$pop.j, 2, sep = "\n")
baSumm$depth.j = factor(baSumm$depth.j)
baSumm$depth.j = factor(baSumm$depth.j, levels = levels(baSumm$depth.j)[c(2, 1)])
#All sites (excluding self retention)
baMeans = baSumm %>% filter(pop.i != pop.j) %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Global")

#mesophotic sources
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Mesophotic") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Mesophotic Source") %>% bind_rows(baMeans, .)

#shallow sources
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Shallow") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Shallow Source") %>% bind_rows(baMeans, .)

#mesophotic sinks
baMeans = baSumm %>% filter(pop.i != pop.j, depth.i == "Mesophotic") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.)))  %>% mutate(dataset = "Mesophotic Sink") %>% bind_rows(baMeans, .)

#shallow sinks
baMeans = baSumm %>% filter(pop.i != pop.j, depth.i == "Shallow") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.)))  %>% mutate(dataset = "Shallow Sink") %>% bind_rows(baMeans, .)

#mesophotic -> shallow
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Mesophotic", depth.i == "Shallow") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Mesophotic -> Shallow") %>% bind_rows(baMeans, .)

#mesophotic -> mesophotic
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Mesophotic", depth.i == "Mesophotic") %>% summarise(mean = mean(mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Mesophotic -> Mesophotic") %>% bind_rows(baMeans, .)

#shallow -> mesophotic
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Shallow", depth.i == "Mesophotic") %>% summarise(mean = mean(.$mean), sd = sd(.$mean), se = sd(.$mean)/sqrt(nrow(.))) %>% mutate(dataset = "Shallow -> Mesophotic") %>% bind_rows(baMeans, .)

#shallow -> shallow
baMeans = baSumm %>% filter(pop.i != pop.j, depth.j == "Shallow", depth.i == "Shallow") %>% summarise(mean = round(mean(.$mean), 3), sd = round(sd(.$mean), 3), se = round(sd(.$mean)/sqrt(nrow(.)), 3)) %>% mutate(dataset = paste("Shallow -> Shallow")) %>% bind_rows(baMeans, .) %>% relocate(dataset, .before = mean) %>% as.data.frame()

baMeans[,c(2:4)] = baMeans[,c(2:4)] %>% round(3)

baMeans
##                    dataset  mean    sd    se
## 1                   Global 0.037 0.047 0.006
## 2        Mesophotic Source 0.055 0.059 0.011
## 3           Shallow Source 0.020 0.020 0.004
## 4          Mesophotic Sink 0.035 0.057 0.011
## 5             Shallow Sink 0.040 0.036 0.007
## 6    Mesophotic -> Shallow 0.049 0.039 0.010
## 7 Mesophotic -> Mesophotic 0.064 0.079 0.023
## 8    Shallow -> Mesophotic 0.014 0.009 0.002
## 9       Shallow -> Shallow 0.028 0.028 0.008
baMeansTabPub = baMeans %>%
  flextable() %>%
  flextable::compose(part = "header", j = "dataset", value = as_paragraph("Dataset")) %>%
  flextable::compose(part = "header", j = "mean", value = as_paragraph(as_i("m"))) %>%
  flextable::compose(part = "header", j = "sd", value = as_paragraph("SD")) %>%
  flextable::compose(part = "header", j = "se", value = as_paragraph("SEM")) %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::fontsize(size = 10, part = "all") %>%
  flextable::bold(part = "header") %>%
  flextable::align(align = "left", part = "all") %>%
  flextable::autofit()

table3 = read_docx()
table3 = body_add_flextable(table3, value = baMeansTabPub)
print(table3, target = "../tables/table3.docx")

baMeansTabPub
baSumm$mean = sprintf('%.3f', baSumm$mean)
baSumm$mean2 = baSumm$mean
baSumm$hpdLow = sprintf('%.3f', baSumm$hpdLow)
baSumm$hpdHigh = sprintf('%.3f', baSumm$hpdHigh)

migrateA = ggplot(data = baSumm, aes(pop.i, pop.j))+
  geom_tile(data = subset(baSumm, subset = baSumm$mean2>0.65), aes(fill = as.numeric(as.character(mean2))), color = "white") +
  geom_segment(data = baSumm, aes(x = 0.4755, xend = -0.625, y = pop.j, yend = pop.j, color = pop.j), size = 14) +
  geom_segment(data = baSumm, aes(x = pop.i, xend = pop.i, y = 0.45, yend = -0.425, color = pop.i), size = 32.25) +
  scale_color_manual(values = flPal[c(1:4, 1:4)], guide = NULL) +
  scale_fill_gradient(low = "grey70", high = "grey5", limit = c(0.65, 0.95), space = "Lab", name = expression(paste(italic("m"))),  guide = "colourbar", breaks = c(0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95)) +
  guides(fill = guide_colorbar(ticks.colour = "black")) +
  new_scale("fill") +
  geom_tile(data = subset(baSumm, subset = baSumm$mean<0.65), aes(fill = as.numeric(as.character(mean))), color = "white") +
  scale_fill_gradientn(colours = paletteer_c("viridis::viridis", n = 10)[3:10], limit = c(0,0.25), space = "Lab", name = expression(paste(italic("m"))), na.value = "transparent",  guide = "colourbar", values = c(0, 0.05, 0.1, 0.15, 0.2,0.5,0.75,1)) +
  geom_text(data = baSumm, aes(x = pop.i, y = pop.j, label = paste(mean, "\n", sep = "")), color = ifelse(baSumm$mean > 0.6, "white", "gray5"), fontface = "bold") +
  geom_text(data = baSumm, aes(x = pop.i, y = pop.j, label = paste("\n(",hpdLow,"–",hpdHigh, ")", sep = "")), color = ifelse(baSumm$mean > 0.6, "white", "gray5"), size = 3.25) +
  labs(x = "Sink", y = "Source") +
  scale_y_discrete(limits = rev(levels(baSumm$pop.i))[c(1:8)], position = "left") +
  coord_cartesian(xlim = c(1, 8), ylim = c(1, 8), clip = "off") +
  theme_minimal()

migrate = migrateA + theme(
  axis.text.x = element_text(vjust = 1, size = 10, hjust = 0.5, color = "gray90"),
  axis.text.y = element_text(size = 10, color = "gray90"),
  axis.title.x = element_text(size = 14),
  axis.title.y = element_text(size = 14),
  panel.grid.major = element_blank(),
  axis.ticks = element_blank(),
  legend.position = "right",
  legend.direction = "vertical",
  legend.title = element_text(size = 14, face = "bold"),
  legend.text = element_text(size = 10)
)

migrate

baSumm$mean = as.numeric(baSumm$mean)
baSumm$hpdLow = as.numeric(baSumm$hpdLow)
baSumm$hpdHigh = as.numeric(baSumm$hpdHigh)

baSummSelf = baSumm %>% filter(pop.i == pop.j) %>% mutate(popdepth = paste(site.i, depth.i)) %>% mutate(retention = mean) %>% dplyr::select(-mean)

fknmsPopsMigrate2 = fknmsSites %>% group_by(site, depthZone, siteID) %>% dplyr::summarise(latDD = first(latDD), longDD = first(longDD)) %>% dplyr::filter(siteID %in% c("Ian's Lumps Site 52", "Site 48", "Site 47", "Site 45", "Site 35/36", "Site 37", "Site 39", "Site 19")) %>% dplyr::select(-site) %>% droplevels() %>% mutate(popdepth = paste(site, depthZone)) %>% as.data.frame() %>% slice(-5) %>% left_join(dplyr::select(baSummSelf, popdepth, retention))
## `summarise()` has grouped output by 'site', 'depthZone'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `site`
## Joining, by = "popdepth"
fknmsPopsMigrate = fknmsPopsMigrate2[c(1,2,4,3,5:8),]

migratePal = c("Upper Keys" = flPal[1], "Lower Keys" = flPal[2], "Tortugas Bank" = flPal[3], "Riley's Hump" = flPal[4])

lines = c("Shallow" = 5, "Mesophotic" = 1)

baMapData = dplyr::select(baSumm, -mean2) %>% left_join(dplyr::select(fknmsPopsMigrate,-retention,-popdepth), by = c("site.i" = "site", "depth.i" = "depthZone")) %>% left_join(dplyr::select(fknmsPopsMigrate,-retention,-popdepth),, by = c("site.j" = "site", "depth.j" = "depthZone"), suffix = c(".i", ".j")) %>% filter(mean >= 0.02)

for(x in 1:nrow(baMapData)) {
  if (baMapData$pop.i[x] == baMapData$pop.j[x]) {
    baMapData$latDD.i[x] = NA;
    baMapData$latDD.j[x] = NA;
    baMapData$longDD.i[x] = NA;
    baMapData$longDD.j[x] = NA;
    baMapData$mean[x] = NA;
    baMapData$median[x] = NA
  }
}

migrateMap = ggplot() +
  geom_sf(data = florida, fill = "white", size = 0.25) +

#SHALLOW SOURCES
#ggtitle("Shallow sources") +
  geom_curve(data = baMapData[2,], aes(x = longDD.j, y = latDD.j, xend = longDD.i, yend = latDD.i-0.03, color = site.j, linetype = depth.j, size = mean), alpha = 0.7, arrow = arrow(type = "open", length = unit(0.03, "npc")), curvature = -3) +
    geom_curve(data = baMapData[7,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 0.25, na.rm = TRUE) +
  geom_curve(data = baMapData[11,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.02, yend = latDD.i-0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 0.2, na.rm = TRUE) +
  geom_curve(data = baMapData[13,], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.04, yend = latDD.i+0.01, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 1.25, na.rm = TRUE) +
  geom_curve(data = baMapData[16,], aes(x = longDD.j, y = latDD.j, xend = longDD.i +0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.25, na.rm = TRUE) +
  geom_curve(data = baMapData[20,], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.04, yend = latDD.i+0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +
  geom_curve(data = baMapData[25,], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.05, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +

#MESO SOURCES
#ggtitle("Mesophotic sources") +
geom_curve(data = baMapData[4,], aes(x = longDD.j, y = latDD.j, xend = longDD.i, yend = latDD.i+0.03, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -3, na.rm = TRUE) +
geom_curve(data = baMapData[c(3, 6),], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.02, yend = latDD.i-0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +
geom_curve(data = baMapData[8,], aes(x = longDD.j, y = latDD.j, xend = longDD.i, yend = latDD.i+0.03, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 0.2, na.rm = TRUE) +
geom_curve(data = baMapData[10,], aes(x = longDD.j, y = latDD.j, xend = longDD.i +0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.25, na.rm = TRUE) +
geom_curve(data = baMapData[12,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.02, yend = latDD.i-0.04, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.35, na.rm = TRUE) +
geom_curve(data = baMapData[15,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.02, yend = latDD.i+0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +
geom_curve(data = baMapData[17,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -.35, na.rm = TRUE) +
geom_curve(data = baMapData[19,], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.04, yend = latDD.i+0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +
geom_curve(data = baMapData[21,], aes(x = longDD.j, y = latDD.j, xend = longDD.i, yend = latDD.i-0.02, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -9, na.rm = TRUE) +
geom_curve(data = baMapData[27,], aes(x = longDD.j, y = latDD.j, xend = longDD.i+0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 2, na.rm = TRUE) +
geom_curve(data = baMapData[c(24, 29),], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.03, yend = latDD.i, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = 0.2, na.rm = TRUE) +
geom_curve(data = baMapData[c(23, 28),], aes(x = longDD.j, y = latDD.j, xend = longDD.i-0.04, yend = latDD.i+0.01, color = site.j, linetype = depth.j, size = mean), arrow = arrow(type = "open", length = unit(0.03, "npc")), alpha = 0.7,  curvature = -0.2, na.rm = TRUE) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = migratePal, guide = NULL) +
  scale_shape_manual(values = c(24, 25), name = "Depth") +
  scale_size(range = c(0.5, 2), breaks = c(0.02,0.06,0.1,0.14,0.18,0.22),name = expression(paste("Migration (", italic("m"), ")", sep = "")), guide = guide_legend(ncol = 3, order = 5)) +
  new_scale("size") +
  geom_point(data = fknmsPopsMigrate, aes(x = longDD, y = latDD, fill = site, shape = depthZone, size = retention)) +
  scale_size(range = c(1.5, 6.5), limit = c(0.68, 0.95), name = expression(paste("Self-recruitment (",italic("m"), ")", sep ="")), breaks = c(0.7, 0.75, 0.8, 0.85, 0.9, 0.95), 
             guide = guide_legend(override.aes = list(shape = 24, fill = "gray80"), ncol = 3, order = 3)) +
  scale_linetype_manual(values = lines, name = "Source depth", guide = guide_legend(order = 4)) +
  coord_sf(xlim = c(-83.1, -80.25), ylim = c(24.3, 25.3)) +
  scale_x_continuous(breaks = c(seq(-84, -80, by = .5))) +
  scale_y_continuous(breaks = c(seq(24, 26, by = .2))) +
  annotation_scale(location = "br") +
  annotation_north_arrow(location = "br", which_north = "true", style = north_arrow_minimal(), pad_x = unit(-0.25, "cm") , pad_y = unit(0.75, "cm")) +
  guides(fill = guide_legend(override.aes = list(shape = 22, color = NA, size = 4),ncol = 2, order = 1), shape = guide_legend(override.aes = list(size = 3), order = 2) ) +
  theme_bw() +
  theme(panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        axis.title = element_blank(),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        plot.title = element_text(color = "black", size = 14),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.box = "horizontal",
        legend.key = element_blank(),
        legend.background = element_blank()
  )

migrateMap

Putting the plots together into a single figure panel

migrationPlots = (migrate/(migrateMap) + theme(legend.box.margin = margin(-20, 0, 0, 0))) + plot_layout(heights = c(1, 1.1)) + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 20))

ggsave("../figures/figure6.png", plot = migrationPlots, width = 27, height = 25, units = "cm", dpi = 300)

ggsave("../figures/figure6.svg", plot = migrationPlots, width = 27, height = 25, units = "cm", dpi = 300)

S. intersepta algal symbiont community structure


Now let’s examine algal symbiont communities with the results of SymPortal analysis of Symbiodiniaceae ITS2 sequences.

ITS2 Data

How many raw reads?

rawItsReads = read.delim("../data/ITS2/sintItsReadCounts", header = FALSE)
colnames(rawItsReads) = c("sample", "reads")

rawItsReads$sample = gsub("_S.*", "", rawItsReads$sample)
rawItsReads = rawItsReads %>% group_by(sample) %>% summarise(reads = first(reads))

head(rawItsReads)
## # A tibble: 6 × 2
##   sample reads
##   <chr>  <int>
## 1 SFK001 34806
## 2 SFK002 38328
## 3 SFK003 18122
## 4 SFK004 52311
## 5 SFK005 30791
## 6 SFK006 40545
#total reads
sum(rawItsReads$reads)
## [1] 7385411
#average reads/sample
(sum(rawItsReads$reads)/nrow(rawItsReads))
## [1] 33723.34
its2Seqs = read.delim("../data/ITS2/148_20210301_DBV_20210401T112728.seqs.absolute.abund_CLEAN.txt", header = TRUE)
its2Profs = read.csv("../data/ITS2/148_20210301_DBV_20210401T112728.profiles.absolute.abund_CLEAN.csv", header = TRUE, check.names = FALSE) 

head(its2Seqs)
##   Sample Symbiodinium  A3 A3b A3at A3ax X43947_A X34778_A X495083_A X36534_A X34149_A
## 1 SFK115            0  12   0    0    0        0        0         0        0        0
## 2 SFK022            0   7   0    0    0        0        0         0        0        0
## 3 SFK025           28 242   7    7    0        6        0         0        7        0
## 4 SFK095            0  14   0    0    0        0        0         0        0        0
## 5 SFK170            0  26   0    0    0        0        0         0        0        0
## 6 SFK175            0   5   0    0    0        0        0         0        0        0
##   X50854_A A3av A3s A3q X33981_A X1402229_A A3au A3aw X50835_A X34175_A X33953_A A3r
## 1        0    0   0   0        0          0    0    0        0        0        0   0
## 2        0    0   0   0        0          0    0    0        0        0        0   0
## 3        0    0  17  23        5          0    0    0        0        0        0   6
## 4        0    0   0   0        0          0    0    0        0        0        0   0
## 5        0    0   0   0        0          0    0    0        0        0        0   0
## 6        0    0   0   0        0          0    0    0        0        0        0   0
##   X1402205_A X364481_A X363583_A A4 X1402230_A X50833_A X50842_A X363143_A X72388_A
## 1          0         0         0  0          0        0        0         0        0
## 2          0         0         0  0          0        0        0         0        0
## 3          0         0         0  0          0        0        0        10        0
## 4          0         0         0  0          0        0        0         0        0
## 5          0         0         0  0          0        0        0         0        0
## 6          0         0         0  7          0        0        0         0        0
##   X363606_A X797686_A X22386_A X529468_A X34696_A X363636_A X1402231_A X34151_A
## 1         0         0        0         0        0         0          0        0
## 2         0         0        0         0        0         0          0        0
## 3         0         0        0         0        0         0          0        0
## 4         0         0        0         0        0         0          0        0
## 5         0         0        0         0        0         0          0        0
## 6         0         0        0         0        0         0          0        0
##   X364459_A X363570_A X363645_A X363578_A X1402232_A X364267_A X363625_A X50845_A
## 1         0         0         0         0          0         0         0        0
## 2         0         0         0         0          0         0         0        0
## 3         0         0         0         0          0         0         0        0
## 4         0         0         0         0          0         0         0        0
## 5         0         0         0         0          0         0         0        0
## 6         0         0         0         0          0         0         0        0
##   X363142_A X1402267_A X22415_A X363639_A X905679_A X363598_A X363706_A X363685_A
## 1         0          0        0         0         0         0         0         0
## 2         0          0        0         0         0         0         0         0
## 3         0          0        0         0         0         0         0         0
## 4         0          0        0         0         0         0         0         0
## 5         0          0        0         0         0         0         0         0
## 6         0          0        0         0         0         0         0         0
##   X1402206_A X1402233_A X1402235_A X43753_A X500385_A X1402234_A A3d X1402202_A
## 1          0          0          0        0         0          0   0          0
## 2          0          0          0        0         0          0   0          0
## 3          0          0          0        0         0          0   0          0
## 4          0          0          0        0         0          0   0          0
## 5          0          0          0        0         0          0   0          0
## 6          0          0          0        0         0          0   0          0
##   X1402236_A X364620_A X363593_A X367833_A X22400_A X50850_A X22463_A X363687_A
## 1          0         0         0         0        0        0        0         0
## 2          0         0         0         0        0        0        0         0
## 3          0         0         0         0        0        0        0         0
## 4          0         0         0         0        0        0        0         0
## 5          0         0         0         0        0        0        0         0
## 6          0         0         0         0        0        0        0         0
##   X693524_A X37988_A X66961_A X22436_A X363590_A X22426_A X45527_A A6b X36953_A
## 1         0        0        0        0         0        0        0   0        0
## 2         0        0        0        0         0        0        0   0        0
## 3         0        0        0        0         0        0        0   0        0
## 4         0        0        0        0         0        0        0   0        0
## 5         0        0        0        0         0        0        0   0        0
## 6         0        0        0        0         0        0        0   0        0
##   X363563_A X37985_A X693526_A X22451_A X33927_A X1402203_A A4.3 X35200_A X22392_A
## 1         0        0         0        0        0          0    0        0        0
## 2         0        0         0        0        0          0    0        0        0
## 3         0        0         0        0        0          0    0        0        0
## 4         0        0         0        0        0          0    0        0        0
## 5         0        0         0        0        0          0    0        0        0
## 6         0        0         0        0        0          0    0        0        0
##   X65140_A X1402245_A X364567_A X364639_A X1402216_A X22464_A X22429_A X49905_A
## 1        0          0         0         0          0        0        0        0
## 2        0          0         0         0          0        0        0        0
## 3        0          0         0         0          0        0        0        0
## 4        0          0         0         0          0        0        0        0
## 5        0          0         0         0          0        0        0        0
## 6        0          0         0         0          0        0        0        0
##   X49571_A X29211_A X36825_A X364218_A X1402237_A X363617_A X73521_A X364172_A
## 1        0        0        0         0          0         0        0         0
## 2        0        0        0         0          0         0        0         0
## 3        0        0        0         0          0         0        0         0
## 4        0        0        0         0          0         0        0         0
## 5        0        0        0         0          0         0        0         0
## 6        0        0        0         0          0         0        0         0
##   X363654_A X1402274_A X49906_A X37990_A X50843_A X363674_A X363646_A X33878_A
## 1         0          0        0        0        0         0         0        0
## 2         0          0        0        0        0         0         0        0
## 3         0          0        0        0        0         0         0        0
## 4         0          0        0        0        0         0         0        0
## 5         0          0        0        0        0         0         0        0
## 6         0          0        0        0        0         0         0        0
##   X1402268_A X1402282_A X22444_A X373280_A X1402238_A X366219_A X69439_A Breviolum B5
## 1          0          0        0         0          0         0        0         0  0
## 2          0          0        0         0          0         0        0         0  0
## 3          0          0        0         0          0         0        0         0  0
## 4          0          0        0         0          0         0        0         0  0
## 5          0          0        0         0          0         0        0         0  0
## 6          0          0        0         0          0         0        0         0  0
##   B18c B18b X1402208_B X1402209_B X45548_B X1402210_B X43411_B X37534_B X1402211_B
## 1    0    0          0          0        0          0        0        0          0
## 2    0    0          0          0        0          0        0        0          0
## 3    0    0          0          0        0          0        0        0          0
## 4    0    0          0          0        0          0        0        0          0
## 5    0    0          0          0        0          0        0        0          0
## 6    0    0          0          0        0          0        0        0          0
##   X1402212_B X1402213_B X1160454_B X1402214_B X1402215_B X37591_B B5ai X1402254_B
## 1          0          0          0          0          0        0    0          0
## 2          0          0          0          0          0        0    0          0
## 3          0          0          0          0          0        0    0          0
## 4          0          0          0          0          0        0    0          0
## 5          0          0          0          0          0        0    0          0
## 6          0          0          0          0          0        0    0          0
##   X71511_B X1402256_B X71527_B X1402255_B X1402257_B X71517_B X71509_B X71508_B
## 1        0          0        0          0          0        0        0        0
## 2        0          0        0          0          0        0        0        0
## 3        0          0        0          0          0        0        0        0
## 4        0          0        0          0          0        0        0        0
## 5        0          0        0          0          0        0        0        0
## 6        0          0        0          0          0        0        0        0
##   X1402258_B X71518_B X71510_B X1402259_B X368876_B X1402260_B X50427_B X1402262_B
## 1          0        0        0          0         0          0        0          0
## 2          0        0        0          0         0          0        0          0
## 3          0        0        0          0         0          0        0          0
## 4          0        0        0          0         0          0        0          0
## 5          0        0        0          0         0          0        0          0
## 6          0        0        0          0         0          0        0          0
##   X71525_B X71523_B X71515_B X368004_B X1402261_B X71526_B X71524_B X1402266_B
## 1        0        0        0         0          0        0        0          0
## 2        0        0        0         0          0        0        0          0
## 3        0        0        0         0          0        0        0          0
## 4        0        0        0         0          0        0        0          0
## 5        0        0        0         0          0        0        0          0
## 6        0        0        0         0          0        0        0          0
##   X1402265_B X1402264_B X1402263_B X43545_B X1402252_B X900132_B X1390171_B X1402253_B
## 1          0          0          0        0          0         0          0          0
## 2          0          0          0        0          0         0          0          0
## 3          0          0          0        0          0         0          0          0
## 4          0          0          0        0          0         0          0          0
## 5          0          0          0        0          0         0          0          0
## 6          0          0          0        0          0         0          0          0
##   B1 X1402276_B X1402280_B X38112_B Cladocopium    C3 C1 C16 C3go C3.10 C42.2 C1dl C3gm
## 1  0          0          0        0         761 10037  0   0    0   123     0    0    0
## 2  0          0          0        0        1273 15896  0   0    0   260     0    0    0
## 3  0          0          0        0         661 11333  0   0    0   115     0    0    0
## 4  0          0          0        0        1019 13663  0   0    0   717     0    0    0
## 5  0          0          0        0        1077 12007  0   0    0   619     0    0    0
## 6  0          0          0        0        1975 21478  0   0    0   385     0    0  110
##   C3gl C3hb C3gr C16b X110271_C X334025_C C3gk C1dk X22330_C X11408_C X18596_C X21897_C
## 1    0    0    0    0         0         0    0    0       58        0        0       90
## 2    0    0    0    0       164       158    0    0      150        0        0      157
## 3    0    0    0    0        67        68    0    0       85        0        0       91
## 4    0    0    0    0         0         0    0    0        0        0        0      227
## 5    0    0    0    0       102       120    0    0       83        0        0       93
## 6    0    0    0    0       282       260    0    0      210        0        0      274
##   C6c C3gq C3gp X65808_C C3gn C15hx C3dw C1cy X1402187_C X20795_C C93.1 X65703_C
## 1   0    0    0       54    0     0    0    0        181        0     0      144
## 2   0    0    0        0    0     0    0    0        120        0     0        0
## 3   0    0    0        0    0     0    0    0          0        0     0       68
## 4   0    0    0       88    0     0    0    0        133        0   229       96
## 5   0    0    0      106    0     0    0    0        131        0    81      120
## 6   0    0    0      166    0     0    0    0        303        0     0        0
##   X385070_C C3ge X1402188_C X1372_C X3238_C X95094_C C3hc X24879_C X91373_C X3699_C
## 1         0    0        107       0       0        0    0        0        0       0
## 2         0    0         93       0     135        0    0      188        0     108
## 3         0    0         76       0       0        0    0      100        0       0
## 4         0    0        120     185       0        0    0        0        0       0
## 5         0    0         89      82       0        0    0        0        0       0
## 6         0    0        128       0       0        0    0      120        0       0
##   C3gt C3dz X20934_C C1af C3gs X25557_C X40208_C X470358_C X40209_C X1402193_C X2239_C
## 1    0    0        0    0    0        0        0         0        0          0       0
## 2    0    0        0    0    0        0        0         0        0          0       0
## 3    0    0        0    0    0        0        0         0        0          0       0
## 4    0    0        0    0    0      170        0         0        0          0       0
## 5    0    0        0    0    0        0        0         0        0          0       0
## 6    0    0      232    0    0        0        0         0        0          0       0
##   X1402195_C C16a X17495_C X17534_C X2097_C X40211_C X93722_C C1v X40207_C X40212_C
## 1          0    0        0        0       0        0        0   0        0        0
## 2          0    0        0        0       0        0        0   0        0        0
## 3         55    0        0        0       0        0        0   0        0        0
## 4          0    0        0        0       0        0        0   0        0        0
## 5          0    0        0        0       0        0        0   0        0        0
## 6          0    0        0        0       0        0        0   0        0        0
##   X1402196_C X1402197_C X9944_C X1402198_C X1402219_C X470998_C X54162_C X22574_C
## 1          0          0       0          0          0         0        0        0
## 2          0          0       0          0          0         0        0        0
## 3          0          0       0          0          0         0        0        0
## 4          0          0       0          0          0         0        0        0
## 5          0          0      94          0          0         0        0        0
## 6          0          0       0          0          0         0        0        0
##   X20921_C X33343_C X1402200_C X25492_C X1402218_C X3240_C X2037_C X85729_C X5371_C
## 1        0        0          0        0          0       0       0        0       0
## 2        0        0          0        0          0       0       0        0       0
## 3        0        0          0       51          0       0       0        0       0
## 4        0        0          0        0          0       0       0        0       0
## 5        0        0          0       57          0       0       0        0       0
## 6        0        0          0        0          0       0       0        0       0
##   X1402225_C X909389_C X1402207_C C3ag X2428_C X1402220_C X4062_C X103828_C X1402199_C
## 1          0         0          0    0       0          0       0       127          0
## 2          0         0          0    0       0          0       0         0          0
## 3          0         0          0    0       0          0       0         0          0
## 4          0         0          0    0       0          0       0         0          0
## 5          0         0          0    0       0          0       0         0          0
## 6          0         0          0    0       0          0       0         0          0
##   X1398518_C X90670_C X1402204_C X1402227_C X1402248_C C6b X1402247_C X1402194_C X870_C
## 1          0        0          0          0          0   0          0          0      0
## 2          0        0          0          0          0   0          0          0      0
## 3          0        0          0          0          0   0          0          0      0
## 4          0        0          0          0          0   0          0          0      0
## 5          0        0          0          0          0   0          0          0      0
## 6          0        0          0          0          0   0          0          0      0
##   X71029_C C3ga X91285_C X1402192_C X1402244_C C1bz X18746_C X1402228_C X694_C C3i
## 1        0    0        0          0          0    0        0          0      0   0
## 2        0    0        0          0          0    0        0          0      0   0
## 3        0    0        0          0          0    0        0          0      0   0
## 4        0    0        0          0          0    0        0          0      0   0
## 5        0    0        0          0          0    0        0          0      0   0
## 6        0    0        0          0          0    0        0          0      0   0
##   X1402250_C X1402243_C C21 X1402281_C C3ca X9108_C X11201_C X11191_C X7821_C
## 1          0          0   0          0    0       0        0        0       0
## 2          0          0   0          0    0       0        0        0       0
## 3          0          0   0          0    0       0        0        0       0
## 4          0          0   0          0    0       0        0        0       0
## 5          0          0   0          0    0       0        0        0       0
## 6          0          0   0          0    0       0        0        0       0
##   X1402273_C C3ck X1402240_C X1402249_C X99988_C X1356_C X69324_C X24193_C C3bb C40f
## 1          0    0          0          0        0       0        0        0    0    0
## 2          0    0          0          0        0       0        0        0    0    0
## 3          0    0          0          0        0       0        0        0    0    0
## 4          0    0          0          0        0       0        0        0    0    0
## 5          0    0          0          0        0       0        0        0    0    0
## 6          0    0          0          0        0       0        0        0    0    0
##   X1401572_C X47282_C X16815_C X5726_C X1402270_C X1402221_C C1ap X1402275_C X21673_C
## 1          0        0        0       0          0          0    0          0        0
## 2          0        0        0       0          0          0    0          0        0
## 3          0        0        0       0          0          0    0          0        0
## 4          0        0        0       0          0          0    0          0        0
## 5          0        0        0       0          0          0    0          0        0
## 6          0        0        0       0          0          0    0          0        0
##   X69758_C C1bt X1402269_C X2943_C C70 X1402271_C X42218_C X1402277_C X9807_C C1ai C3t
## 1        0    0          0       0   0          0        0          0       0    0   0
## 2        0    0          0       0   0          0        0          0       0    0   0
## 3        0    0          0       0   0          0        0          0       0    0   0
## 4        0    0          0       0   0          0        0          0       0    0   0
## 5        0    0          0       0   0          0        0          0       0    0   0
## 6        0    0          0       0   0          0        0          0       0    0   0
##   X54249_C X26258_C X1402278_C X1402272_C C1x X40218_C X21804_C X1402279_C C3de
## 1        0        0          0          0   0        0        0          0    0
## 2        0        0          0          0   0        0        0          0    0
## 3        0        0          0          0   0        0        0          0    0
## 4        0        0          0          0   0        0        0          0    0
## 5        0        0          0          0   0        0        0          0    0
## 6        0        0          0          0   0        0        0          0    0
##   X13929_C X23354_C X1402223_C X99010_C X983542_C X3366_C X1402226_C X1402222_C
## 1        0        0          0        0         0       0          0          0
## 2        0        0          0        0         0       0          0          0
## 3        0        0          0        0         0       0          0          0
## 4        0        0          0        0         0       0          0          0
## 5        0        0          0        0         0       0          0          0
## 6        0        0          0        0         0       0          0          0
##   X42529_C X2152_C X62532_C X4558_C X2427_C X1829_C C3fo X18793_C X11809_C X31248_C
## 1        0       0        0       0       0       0    0        0        0        0
## 2        0       0        0       0       0       0    0        0        0        0
## 3        0       0        0       0       0       0    0        0        0        0
## 4        0       0        0       0       0       0    0        0        0        0
## 5        0       0        0       0       0       0    0        0        0        0
## 6        0       0        0       0       0       0    0        0        0        0
##   X1402241_C X816_C X921460_C C3ao C3an C3cn X3241_C X103581_C X21093_C X1402224_C
## 1          0      0         0    0    0    0       0         0        0          0
## 2          0      0         0    0    0    0       0         0        0          0
## 3          0      0         0    0    0    0       0         0        0          0
## 4          0      0         0    0    0    0       0         0        0          0
## 5          0      0         0    0    0    0       0         0        0          0
## 6          0      0         0    0    0    0       0         0        0          0
##   X18813_C X1402201_C X22869_C X23865_C C6a C1j X1402239_C X3434_C X22178_C X17016_C
## 1        0          0        0        0   0   0          0       0        0        0
## 2        0          0        0        0   0   0          0       0        0        0
## 3        0          0        0        0   0   0          0       0        0        0
## 4        0          0        0        0   0   0          0       0        0        0
## 5        0          0        0        0   0   0          0       0        0        0
## 6        0          0        0        0   0   0          0       0        0        0
##   X1402242_C X42518_C X54160_C X873_C C50f X1402246_C X1390080_C X113247_C X99987_C
## 1          0        0        0      0    0          0          0         0        0
## 2          0        0        0      0    0          0          0         0        0
## 3          0        0        0      0    0          0          0         0        0
## 4          0        0        0      0    0          0          0         0        0
## 5          0        0        0      0    0          0          0         0        0
## 6          0        0        0      0    0          0          0         0        0
##   X3601_C X1866_C X1402251_C X2895_C X9153_C C3fn X866_C X864_C X18159_C X990_C
## 1       0       0          0       0       0    0      0      0        0      0
## 2       0       0          0       0       0    0      0      0        0      0
## 3       0       0          0       0       0    0      0      0        0      0
## 4       0       0          0       0       0    0      0      0        0      0
## 5       0       0          0       0       0    0      0      0        0      0
## 6       0       0          0       0       0    0      0      0        0      0
##   X1402189_C X21205_C C3fc X871_C X1402190_C X8117_C X55844_C X54218_C X51874_C
## 1          0        0    0      0          0       0        0        0        0
## 2          0        0    0      0          0       0        0        0        0
## 3          0        0    0      0          0       0        0        0        0
## 4          0        0    0      0          0       0        0        0        0
## 5          0        0    0      0          0       0        0        0        0
## 6          0        0    0      0          0       0        0        0        0
##   X874099_C X27927_C C65b X46049_C X37410_C X28411_C D1 G3l G3b
## 1         0        0    0        0        0        0  0   0   0
## 2         0        0    0        0        0        0  0   0   0
## 3         0        0    0        0        0        0  0   0   0
## 4         0        0    0        0        0        0  0   0   0
## 5         0        0    0        0        0        0  0   0   0
## 6         0        0    0        0        0        0  0   0   0
sum(its2Seqs[,c(2:ncol(its2Seqs))])
## [1] 3280586
sum(its2Profs[,c(2:ncol(its2Profs))])
## [1] 2728837
its2SeqsGen = its2Seqs %>% rowwise() %>%  summarise(sample = Sample, symbiodinium = sum(c_across(2:121)), breviolum = sum(c_across(122:176)), cladocopium = sum(c_across(177:405)), durusdinium = sum(c_across(406)), gerakladium = sum(c_across(406:407)))

round(sum(its2SeqsGen$symbiodinium)/sum(its2SeqsGen[,-1])*100, 2)
## [1] 19.67
round(sum(its2SeqsGen$breviolum)/sum(its2SeqsGen[,-1])*100, 2)
## [1] 0.29
round(sum(its2SeqsGen$cladocopium)/sum(its2SeqsGen[,-1])*100, 2)
## [1] 80.03
round(sum(its2SeqsGen$durusdinium)/sum(its2SeqsGen[,-1])*100, 4)
## [1] 0.0002
round(sum(its2SeqsGen$gerakladium)/sum(its2SeqsGen[,-1])*100, 4) 
## [1] 0.002
its2ProfsGen = its2Profs %>% rowwise() %>%  summarise(sample = Sample, symbiodinium = sum(c_across(2:7)), breviolum = sum(c_across(8:10)), cladocopium = sum(c_across(11:20)))

round(sum(its2ProfsGen$symbiodinium)/sum(its2ProfsGen[,-1])*100, 2)
## [1] 21.33
round(sum(its2ProfsGen$breviolum)/sum(its2ProfsGen[,-1])*100, 2)
## [1] 0.22
round(sum(its2ProfsGen$cladocopium)/sum(its2ProfsGen[,-1])*100, 2)
## [1] 78.45

Read in SymPortal outputs for ITS2 type profiles

stephanocoeniaMetaData = read.csv("../data/stephanocoeniaMetaData.csv", header = TRUE, check.names = FALSE)[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select(c(Sample = tubeID, site, depthM, depthZone))

its2Profs = read.csv("../data/ITS2/148_20210301_DBV_20210401T112728.profiles.absolute.abund_CLEAN.csv", header = TRUE, check.names = FALSE) 

its2Profs = stephanocoeniaMetaData %>% right_join(its2Profs) %>% arrange(Sample) 
## Joining, by = "Sample"
its2Profs$site = factor(its2Profs$site)
its2Profs$site = factor(its2Profs$site, levels(its2Profs$site)[c( 2, 3, 1, 4)])
its2Profs$depthZone = factor(its2Profs$depthZone)
its2Profs$depthZone = factor(its2Profs$depthZone, levels(its2Profs$depthZone)[c(2, 1)])

its2Profs = its2Profs %>% arrange(site, depthZone,   desc(`C3/C3.10`),  desc(`C1/C3-C42.2-C1dl-C3gl-C3gm-C3gk`), desc(`C3-C1-C3.10`), desc(`C3-C1dk-C15hx`), desc(`C3-C3go-C6c-C3gq-C3gp-C3gn-C3dw`), desc(`C16/C3-C16b`), desc(`C3-C3hb-C3ge-C3hc-C1dk`), desc(`C3-C3gr-C3gt-C3gs-C3.10`), desc(`C3/C1`),desc(`A3-A3b-A3at-A3ax`), desc(`A3-A3at-A3b-A3q-A3s`), desc(`A3-A3s-A3q`), desc(`A3`), desc(`A3-A3b-A3av-A3au-A3aw`),desc(`A4`), desc(`C3`), desc(`B18b`), desc(`B18c`), desc(`B5`))

sampleCounts = plyr::count(its2Profs, c('site','depthZone'))
meltedList = reshape2::melt(lapply(sampleCounts$freq,function(x){c(1:x)}))
its2Profs$barPlotOrder = meltedList$value
its2Profs = its2Profs[c(1,ncol(its2Profs),2:(ncol(its2Profs)-1))]

head(its2Profs)



Preparing ITS2 type profiles for plotting

its2ProfsPerc = its2Profs
its2ProfsPerc$sum = apply(its2ProfsPerc[, c(6:length(its2ProfsPerc[1,]))], 1, function(x) {
sum(x, na.rm = T)
})

its2ProfsPerc = cbind(its2ProfsPerc[, c(1:5)], (its2ProfsPerc[, c(6:(ncol(its2ProfsPerc)-1))]
/ its2ProfsPerc$sum))
head(its2ProfsPerc)
##   Sample barPlotOrder         site depthM depthZone A3-A3b-A3at-A3ax
## 1 SFK068            1 Riley's Hump   27.4   Shallow                0
## 2 SFK091            2 Riley's Hump   26.2   Shallow                0
## 3 SFK073            3 Riley's Hump   26.2   Shallow                0
## 4 SFK084            4 Riley's Hump   26.2   Shallow                0
## 5 SFK083            5 Riley's Hump   26.2   Shallow                0
## 6 SFK082            6 Riley's Hump   26.5   Shallow                0
##   A3-A3at-A3b-A3q-A3s A3-A3s-A3q A3 A3-A3b-A3av-A3au-A3aw A4 B18b B18c B5 C3/C3.10
## 1                   0          0  0                     0  0    0    0  0        0
## 2                   0          0  0                     0  0    0    0  0        0
## 3                   0          0  0                     0  0    0    0  0        0
## 4                   0          0  0                     0  0    0    0  0        0
## 5                   0          0  0                     0  0    0    0  0        0
## 6                   0          0  0                     0  0    0    0  0        0
##   C1/C3-C42.2-C1dl-C3gl-C3gm-C3gk C3-C1-C3.10 C3-C1dk-C15hx
## 1                               1           0             0
## 2                               1           0             0
## 3                               1           0             0
## 4                               1           0             0
## 5                               1           0             0
## 6                               1           0             0
##   C3-C3go-C6c-C3gq-C3gp-C3gn-C3dw C16/C3-C16b C3-C3hb-C3ge-C3hc-C1dk
## 1                               0           0                      0
## 2                               0           0                      0
## 3                               0           0                      0
## 4                               0           0                      0
## 5                               0           0                      0
## 6                               0           0                      0
##   C3-C3gr-C3gt-C3gs-C3.10 C3/C1 C3
## 1                       0     0  0
## 2                       0     0  0
## 3                       0     0  0
## 4                       0     0  0
## 5                       0     0  0
## 6                       0     0  0
# check that all proportions add up to 1
apply(its2ProfsPerc[, c(6:(ncol(its2ProfsPerc)))], 1, function(x) {
sum(x, na.rm = T)
})
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [42] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [124] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [206] 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Everything looks good and is ready to plot

gssProf = otuStack(its2ProfsPerc, count.columns = c(6:length(its2ProfsPerc[1, ])),
 condition.columns = c(1:5)) %>% filter(otu != "summ") %>% droplevels() # remove summ rows

levels(gssProf$otu)
levels(gssProf$depthZone)
levels(gssProf$site)



Consruct ITS2 type profile barplot

getPaletteA = colorRampPalette(c("#006994","#FFFFFF"))
getPaletteB = colorRampPalette(c("#FFBF46","#FFFFFF"))
getPaletteC = colorRampPalette(c("#085F63","#FFFFFF"))


profPal = c(getPaletteA(8)[1:6], getPaletteB(5)[1:3], getPaletteC(11)[1:10])

popAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5), x2 = c(30.5, 30.5, 30.5, 30.5),
                     y1 = -0.14, y2 = -0.14,
                     site = c("Riley's Hump", "Tortugas Bank", "Lower Keys", "Upper Keys"))

popAnno$site = factor(popAnno$site)
popAnno$site = factor(popAnno$site, levels = levels(popAnno$site)[c(2, 3, 1, 4)])
 
gssProfPlot = gssProf %>% left_join(popAnno, by = "site")

its2ProfsPlotA = ggplot(gssProfPlot, aes(x = barPlotOrder, y = count, fill = otu)) +
  geom_bar(position = "stack", stat = "identity", color = "gray25", size = 0.25) + 
  scale_fill_manual(values = profPal) +
  scale_color_manual(values = rev(flPal)) +
  geom_segment(aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 8) +
  labs(title = expression(italic("SymPortal")), fill = expression(paste(italic("ITS2"), " type profile"))) +
  guides(color = "none", fill = guide_legend(ncol = 3, reverse = FALSE)) +
  facet_grid(factor(depthZone) ~ site, scales = "free", switch = "both", space = "free") + # faceting plots by Depth and Site
  scale_x_discrete(expand = c(0.03, 0.03)) +
  scale_y_continuous(expand = c(0.002, 0.002)) +
  coord_cartesian(ylim = c(-.01, 1.01), clip = "off") +
theme_bw()

its2ProfsPlot = its2ProfsPlotA +
theme(plot.title = element_text(),
  panel.grid = element_blank(),
  panel.background = element_rect(fill = "gray25", colour = "gray25"),
  panel.border = element_rect(fill = NA, color = "black", size = 1, linetype = "solid"),
  panel.spacing.x = grid:::unit(0.05, "lines"),
  panel.spacing.y = grid:::unit(0.05, "lines"),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_blank(),
  axis.ticks.y = element_blank(),
  axis.title = element_blank(),
  strip.background.x = element_blank(),
  strip.background.y = element_blank(),
  strip.text = element_text(size = 12),
  strip.text.y.left = element_text(angle = 90),
  strip.text.x.bottom = element_text(vjust = 0, color = "grey90"),
  legend.key.size = unit(0.75, "line"),
  legend.key = element_blank(),
  legend.title = element_text(size = 10, angle = 90),
  legend.text = element_text(size = 8),
  legend.position = "bottom")

its2ProfsPlot



SNP vs ITS2 genera

Pulling genera of Symbiodiniaceae from SNPS and comparing to genera of ITS2 profiles from SymPortal

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("Sample" = tubeID, "site" = site, "depth" = depthZone)

zoox = read.delim("../data/snps/symbionts/zooxReads", header = FALSE, check.names = FALSE)

head(zoox)
##                     V1  V2
## 1 fk_S001.trim.bt2.bam  NA
## 2                 chr1 103
## 3                 chr2  95
## 4                 chr3 135
## 5                 chr4 123
## 6                 chr5  11
#Reconstruct read mapping output into dataframe usable for analysis
zoox$V2[is.na(zoox$V2)] <- as.character(zoox$V1[is.na(zoox$V2)])
zoox$V1 = gsub(pattern = "fk_*", "chr", zoox$V1)
zoox$V2 = gsub(".trim.*", "", zoox$V2)
zoox = zoox %>% filter(zoox$V1 != "*")
zooxLst = split(zoox$V2, as.integer(gl(length(zoox$V2), 20, length(zoox$V2))))

zooxMaps = NULL

for(i in zooxLst){
  zooxMaps = rbind(zooxMaps, data.frame(t(i)))
}

#remove tech reps
zooxMaps = zooxMaps[-c(66, 68, 164, 166, 209, 211),]

#rename columns and samples to match other ITS2 dataframe
zooxMaps$X1 = gsub("fk_S", "SFK", zooxMaps$X1)
zooxMaps$X1 = gsub("\\.[1-3]", "", zooxMaps$X1)
colnames(zooxMaps) = c("Sample",zoox$V1[c(2:20)])

#convert characters to numeric
str(zooxMaps)
## 'data.frame':    220 obs. of  20 variables:
##  $ Sample: chr  "SFK001" "SFK002" "SFK003" "SFK004" ...
##  $ chr1  : chr  "103" "86" "57" "970" ...
##  $ chr2  : chr  "95" "138" "63" "1096" ...
##  $ chr3  : chr  "135" "155" "78" "1471" ...
##  $ chr4  : chr  "123" "182" "34" "1394" ...
##  $ chr5  : chr  "11" "15" "5" "26" ...
##  $ chr6  : chr  "75" "72" "110" "42" ...
##  $ chr7  : chr  "93" "134" "123" "43" ...
##  $ chr8  : chr  "148" "230" "125" "100" ...
##  $ chr9  : chr  "98" "90" "78" "39" ...
##  $ chr10 : chr  "3131" "3973" "3664" "2855" ...
##  $ chr11 : chr  "4468" "5327" "5232" "3649" ...
##  $ chr12 : chr  "4696" "5532" "5511" "3902" ...
##  $ chr13 : chr  "4360" "5070" "5102" "3546" ...
##  $ chr14 : chr  "3750" "4587" "4354" "3294" ...
##  $ chr15 : chr  "549" "590" "554" "458" ...
##  $ chr16 : chr  "35" "121" "24" "87" ...
##  $ chr17 : chr  "34" "70" "25" "38" ...
##  $ chr18 : chr  "25" "61" "21" "34" ...
##  $ chr19 : chr  "3" "6" "1" "11" ...
for(i in c(2:20)){
  zooxMaps[,i] = as.numeric(zooxMaps[,i])
  }

str(zooxMaps)
## 'data.frame':    220 obs. of  20 variables:
##  $ Sample: chr  "SFK001" "SFK002" "SFK003" "SFK004" ...
##  $ chr1  : num  103 86 57 970 73 3600 31 52 91 46 ...
##  $ chr2  : num  95 138 63 1096 89 ...
##  $ chr3  : num  135 155 78 1471 121 ...
##  $ chr4  : num  123 182 34 1394 79 ...
##  $ chr5  : num  11 15 5 26 7 118 15 14 12 14 ...
##  $ chr6  : num  75 72 110 42 52 72 58 78 99 76 ...
##  $ chr7  : num  93 134 123 43 76 76 52 143 166 126 ...
##  $ chr8  : num  148 230 125 100 118 124 99 185 215 191 ...
##  $ chr9  : num  98 90 78 39 50 58 37 107 105 73 ...
##  $ chr10 : num  3131 3973 3664 2855 1518 ...
##  $ chr11 : num  4468 5327 5232 3649 1877 ...
##  $ chr12 : num  4696 5532 5511 3902 2097 ...
##  $ chr13 : num  4360 5070 5102 3546 1802 ...
##  $ chr14 : num  3750 4587 4354 3294 1783 ...
##  $ chr15 : num  549 590 554 458 255 200 500 823 589 865 ...
##  $ chr16 : num  35 121 24 87 79 99 26 27 78 36 ...
##  $ chr17 : num  34 70 25 38 42 57 30 51 58 50 ...
##  $ chr18 : num  25 61 21 34 54 59 17 24 44 39 ...
##  $ chr19 : num  3 6 1 11 2 16 3 0 0 18 ...
#collapse fake chromosomes into representativ genera
zooxMaps$Symbiodinium = rowSums(zooxMaps[2:6])
zooxMaps$Breviolum = rowSums(zooxMaps[7:10])
zooxMaps$Cladcopium = rowSums(zooxMaps[11:16])
zooxMaps$Durusdinium = rowSums(zooxMaps[17:20])

#keep genera totals and turn into proportions for barplot
zooxMaps = zooxMaps[,c(1, 21:24)]
zooxProp = zooxMaps
zooxProp$sum = apply(zooxProp[, c(2:length(zooxProp[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
zooxProp = cbind(zooxProp$Sample, (zooxProp[, c(2:(ncol(zooxProp)-1))]
/ zooxProp$sum))

colnames(zooxProp)[1] = "Sample"

head(zooxProp)
##   Sample Symbiodinium   Breviolum Cladcopium Durusdinium
## 1 SFK001  0.021293088 0.018876527  0.9554076 0.004422761
## 2 SFK002  0.021785998 0.019894852  0.9485608 0.009758312
## 3 SFK003  0.009419339 0.017328405  0.9704304 0.002821827
## 4 SFK004  0.215007591 0.009715897  0.7679028 0.007373672
## 5 SFK005  0.036268921 0.029093768  0.9172400 0.017397287
## 6 SFK006  0.704310476 0.012543713  0.2743652 0.008780599
#Check that all samples total to 1
apply(zooxProp[, c(2:(ncol(zooxProp)))], 1, function(x) {
sum(x, na.rm = T)
})
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  67 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 157 158 159 160 161 162 163 165 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 203 204 205 206 207 208 210 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
#add sample metadata to proportions
snpSym = popData %>% left_join(zooxProp)
## Joining, by = "Sample"

Combining SNP and ITS2 data for comparison of Symbiodiniaceae genera This will allow us to plot individuals in the same order across methods

#sum profiles into genera
symGenera = its2Profs
symGenera$itsSymbiodinium = rowSums(symGenera[6:11])
symGenera$itsBreviolum = rowSums(symGenera[12:14])
symGenera$itsCladcopium = rowSums(symGenera[15:24])
symGenera$itsDurusdinium = 0  

symGenera = symGenera %>% dplyr::select(Sample, barPlotOrder, itsSymbiodinium, itsBreviolum, itsCladcopium, itsDurusdinium)

#convert to proportions
symGenera$sum = apply(symGenera[, c(3:length(symGenera[1,]))], 1, function(x) {
sum(x, na.rm = T)
})

symGeneraProp = cbind(symGenera$Sample, symGenera[, c(3:(ncol(symGenera)-1))]
/ symGenera$sum)

colnames(symGeneraProp)[1] = "Sample"

#Check that all samples total to 1
apply(symGeneraProp[,c(2:5)], 1, function(x) {
sum(x, na.rm = T)
})
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [42] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [124] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [206] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#construct combined dataframe
symGenera = symGenera %>% dplyr::select(Sample, barPlotOrder) %>% left_join(snpSym) %>%  left_join(symGeneraProp) 
## Joining, by = "Sample"
## Joining, by = "Sample"
symGenera$depth = factor(symGenera$depth)
symGenera$depth = factor(symGenera$depth, levels = levels(symGenera$depth)[c(2, 1)])

symGenera$site = factor(symGenera$site)
symGenera$site = factor(symGenera$site, levels = levels(symGenera$site)[c(2, 3, 1, 4)])

# symGenera$barPlotOrder[is.na(symGenera$barPlotOrder)] <- 30

#turn into melted dataframe with otustack() and remove "summ" rows
gssSym = otuStack(symGenera, count.columns = c(5:length(symGenera[1, ])),
 condition.columns = c(1:4)) %>% filter(otu != "summ") %>% droplevels()

#check that levels are correct/ordered
levels(gssSym$otu)
## [1] "Symbiodinium"    "Breviolum"       "Cladcopium"      "Durusdinium"    
## [5] "itsSymbiodinium" "itsBreviolum"    "itsCladcopium"   "itsDurusdinium"
levels(gssSym$depth)
## [1] "Shallow"    "Mesophotic"
levels(gssSym$site)
## [1] "Riley's Hump"  "Tortugas Bank" "Lower Keys"    "Upper Keys"


Creating Symbiodiniaceae genera relative proportion barplots

SNPs:

colPalZoox = c("#247EA3", "#FFBF46", "#6A9FA1", "Purple3")

popAnno = data.frame(x1 = c(0.25, 0.25, 0.25, 0.25), x2 = c(30.75, 30.75, 30.75, 30.75),
                     y1 = -0.15, y2 = -0.15, site = c("Riley's Hump", "Tortugas Bank", "Lower Keys", "Upper Keys"))

popAnno$site = factor(popAnno$site)
popAnno$site = factor(popAnno$site, levels = levels(popAnno$site)[c(2, 3, 1, 4)])

gssSymPlot = gssSym %>% left_join(popAnno, by = "site")

zooxSNPA = ggplot(data = subset(gssSymPlot, subset = otu %in% c("Symbiodinium", "Breviolum", "Cladcopium", "Durusdinium" )), aes(x = barPlotOrder, y = count, fill = otu, order = barPlotOrder)) +
  geom_point(aes(x=0, y=0, fill = otu), shape = 22, size = 0) +
  geom_bar(stat = "identity", position = "stack", colour = "grey25", width = 1, size = 0.2, show.legend = FALSE) +
  xlab("Population") +
  scale_x_discrete(expand = c(0.035, 0.035)) +
  scale_y_continuous(expand = c(0.002, 0.002)) +
  scale_color_manual(values = rev(flPal)) +
  geom_segment(aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 7) +
  scale_fill_manual(values = colPalZoox, name = "Symbiodiniaceae genus") +
  coord_cartesian(ylim = c(-.01,1.01), clip = "off") +
  facet_grid(factor(depth) ~ site, drop = TRUE, scales = "free", switch = "both", space = "free") +
  ggtitle("2bRAD") +
  guides(fill = guide_legend(override.aes = list(size = 5), ncol = 1), color = "none")  +
  theme_bw()

zooxSNP = zooxSNPA + theme(plot.title = element_text(),
  panel.grid = element_blank(),
  panel.background = element_rect(fill = "gray25", colour = "grey25"),
  panel.border = element_rect(fill = NA, color = "black", size = 1, linetype = "solid"),
  panel.spacing.x = grid:::unit(0.05, "lines"),
  panel.spacing.y = grid:::unit(0.05, "lines"),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_blank(),
  axis.ticks.y = element_blank(),
  axis.title = element_blank(),
  strip.background.x = element_blank(),
  strip.background.y = element_blank(),
  strip.text = element_text(size = 12),
  strip.text.y.left = element_text(size = 12, angle = 90),
  strip.text.x.bottom = element_text(vjust = -.1, color = "grey90"),
  legend.title = element_text(size = 10),
  legend.text = element_text(size = 8, face = "italic"),
  legend.key.size = unit(0.5, "line"),
  legend.key = element_blank(),
  legend.position = "bottom",
  legend.direction = "vertical",
  legend.box = "horizontal")

# zooxSNP


ITS2:

zooxITSA = ggplot(data = subset(gssSymPlot, subset = !(otu %in% c("Symbiodinium", "Breviolum", "Cladcopium", "Durusdinium" ))), aes(x = barPlotOrder, y = count, fill = otu, order = barPlotOrder)) +
  geom_point(aes(x=0, y=0, fill = otu), shape = 22, size = 0) +
  geom_bar(stat = "identity", position = "stack", colour = "grey25", width = 1, size = 0.2, show.legend = FALSE) +
  xlab("Population") +
  scale_fill_manual(values = colPalZoox, name = "Symbiodiniaceae genus", labels = c("Symbiodinium", "Breviolum", "Cladcopium", "Durusdinium")) +
  scale_color_manual(values = rev(flPal)) +
  geom_segment(aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 7) +
  coord_cartesian(ylim = c(-.01,1.01), clip = "off") +
  scale_x_discrete(expand = c(0.035, 0.035)) +
  scale_y_continuous(expand = c(0.002, 0.002)) +
  facet_grid(factor(depth) ~ site, drop = TRUE, scales = "free", switch = "both", space = "free") +
  guides(fill = guide_legend(override.aes = list(size = 5), ncol = 1), color = "none")  +
  labs(title = expression(italic("ITS2"))) +
  theme_bw()

zooxITS = zooxITSA + theme(plot.title = element_text(),
  panel.grid = element_blank(),
  panel.background = element_rect(fill = "gray25", colour = "grey25"),
  panel.border = element_rect(fill = NA, color = "black", size = 1, linetype = "solid"),
  panel.spacing.x = grid:::unit(0.05, "lines"),
  panel.spacing.y = grid:::unit(0.05, "lines"),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_blank(),
  axis.ticks.y = element_blank(),
  axis.title = element_blank(),
  strip.background.x = element_blank(),
  strip.background.y = element_blank(),
  strip.text = element_text(size = 12),
  strip.text.y.left = element_text(size = 12, angle = 90),
  strip.text.x.bottom = element_text(vjust = -.1, color = "grey90"),
  legend.title = element_text(size = 10),
  legend.text = element_text(size = 8, face = "italic"),
  legend.key = element_blank(),
  legend.key.size = unit(0.1, "line"),
  legend.position = "bottom",
  legend.direction = "vertical",
  legend.box = "horizontal")
 
# zooxITS

Symbiodiniaceae barplots

its2Plots = (its2ProfsPlot + theme(legend.position = "right",  legend.title = element_text(angle = 0)) & guides(color = "none", fill = guide_legend(ncol = 1, reverse = FALSE)))/zooxSNP/zooxITS +
  plot_annotation(tag_levels = "A") &
  theme(plot.tag = element_text(size = 16), 
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "right",
        legend.title = element_text(color = "black", size = 10),
        legend.text = element_text(color = "black", size = 8))


ggsave("../figures/figure7.png", plot = its2Plots, height = 26, width = 20, units = "cm", dpi = 300)
ggsave("../figures/figure7.svg", plot = its2Plots, height = 26, width = 20, units = "cm", dpi = 300)

Comparing the two outputs with procrustes analysis

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("Sample" = tubeID, "site" = site, "depth" = depthZone)

symSnpDf = zooxMaps %>% left_join(popData) %>% relocate(c(site, depth), .after = Sample) %>% filter(!row_number()==131) %>% mutate(dataSet = "SNPs") %>% relocate(dataSet, .after = Sample)
## Joining, by = "Sample"
rownames(symSnpDf) = symSnpDf$Sample

symITS2 = its2Profs
symITS2$Symbiodinium = rowSums(symITS2[6:11])
symITS2$Breviolum = rowSums(symITS2[12:14])
symITS2$Cladcopium = rowSums(symITS2[15:24])
symITS2$Durusdinium = 0  

symITS2Df = symITS2 %>% dplyr::select(Sample, Symbiodinium, Breviolum, Cladcopium, Durusdinium) %>% left_join(popData) %>% relocate(c(site, depth), .after = Sample) %>% arrange(Sample) %>% mutate(dataSet = "ITS2") %>% relocate(dataSet, .after = Sample)
## Joining, by = "Sample"
rownames(symITS2Df) = symITS2Df$Sample

#create distance matrices
symSnpdist = vegdist(decostand(symSnpDf[c(5:ncol(symSnpDf))], method = "normalize"), method = "bray")
    
symITS2dist = vegdist(decostand(symITS2Df[c(5:ncol(symITS2Df))], method = "normalize"), method = "bray")

snpPcOA = cmdscale(symSnpdist, eig = TRUE, x.ret = TRUE)
its2PcOA = cmdscale(symITS2dist, eig = TRUE, x.ret = TRUE)

#procrustes analysis
its2GeneraProcrustes = protest(Y = its2PcOA, X = snpPcOA, choices = c(1, 2), 
permutations = 9999, symmetric = FALSE)

its2GeneraProcrustes
## 
## Call:
## protest(X = snpPcOA, Y = its2PcOA, permutations = 9999, choices = c(1,      2), symmetric = FALSE) 
## 
## Procrustes Sum of Squares (m12 squared):        0.1582 
## Correlation in a symmetric Procrustes rotation: 0.9175 
## Significance:  0.0001 
## 
## Permutation: free
## Number of permutations: 9999
plot(its2GeneraProcrustes, kind = 1)

plot(its2GeneraProcrustes, kind = 2)

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "site" = site, "depth" = depthZone)

symGenProcPlot = procrustes(X = snpPcOA, Y = its2PcOA, choices = c(1, 2), symmetric = FALSE)

symGenProcPlotData = cbind(symGenProcPlot$X, symGenProcPlot$Yrot) %>% as.data.frame()
rownames(symGenProcPlotData) = rownames(symGenProcPlot$X)
colnames(symGenProcPlotData) = c("x1", "y1", "x2", "y2")
symGenProcPlotData$sample = row.names(symGenProcPlotData)
symGenProcPlotData$sample = gsub(pattern = "\\.2", "", symGenProcPlotData$sample)
symGenProcPlotData = symGenProcPlotData %>% left_join(popData) %>% relocate(sample, .before = x1)
## Joining, by = "sample"
symGenProcPlotData$depth = factor(symGenProcPlotData$depth)
symGenProcPlotData$depth = factor(symGenProcPlotData$depth, levels(symGenProcPlotData$depth)[c(2,1)])
symGenProcPlotData$site = factor(symGenProcPlotData$site)
symGenProcPlotData$site = factor(symGenProcPlotData$site, levels(symGenProcPlotData$site)[c(4, 1, 3, 2)])
  
head(symGenProcPlotData)
##   sample          x1           y1         x2            y2          site      depth
## 1 SFK001 -0.11319743  0.009463110 -0.1427982 -0.0009174297 Tortugas Bank Mesophotic
## 2 SFK002 -0.11052338  0.010627368 -0.1264448  0.0067320405 Tortugas Bank Mesophotic
## 3 SFK003 -0.12633478  0.001067354 -0.1427982 -0.0009174297 Tortugas Bank Mesophotic
## 4 SFK004  0.06580863 -0.060023821  0.1332136 -0.0332631646 Tortugas Bank Mesophotic
## 5 SFK005 -0.09426666  0.029200059 -0.1427982 -0.0009174297 Tortugas Bank Mesophotic
## 6 SFK006  0.55445199 -0.083081992  0.4276686 -0.0594228026 Tortugas Bank Mesophotic
#Calculate the angle of rotation, and then the slope of the rotated axis
theta = acos(symGenProcPlot$rotation[1,1]) 
slope = tan(theta)

#Calculate the y-intercepts for rotated axes
symGenProcInt = symGenProcPlot$translation[2] - (slope*symGenProcPlot$translation[1])
symGenProcInt2 = symGenProcPlot$translation[2] - (-1/slope*symGenProcPlot$translation[1])

sintSymGenProcPlotA = ggplot() +
  geom_hline(yintercept = 0, color = "gray90", linetype = 1) +
  geom_vline(xintercept = 0, color = "gray90", linetype = 1) +
  geom_abline(intercept = symGenProcInt, slope = slope, color = "gray75", linetype = 2) +
  geom_abline(intercept = symGenProcInt2, slope = -(1/slope), color = "gray75", linetype = 2) +
    geom_segment(data = symGenProcPlotData, aes(x = x2*(1-symGenProcPlot$scale), y = y2*(1-symGenProcPlot$scale), xend = x1*(1-symGenProcPlot$scale), yend = y1*(1-symGenProcPlot$scale), color = site), alpha = 0.5) +
  geom_point(data = symGenProcPlotData, aes(x = x2*(1-symGenProcPlot$scale), y = y2*(1-symGenProcPlot$scale), fill = site, shape = depth), alpha = 0.5)+
  geom_point(data = symGenProcPlotData, aes(x = x1*(1-symGenProcPlot$scale), y = y1*(1-symGenProcPlot$scale), fill = site, shape = depth), size = 2) +
  annotate(geom = "label", x = 0.172, y = 0.172, label = "           ", size = 10) +
  annotate(geom = "text", x = 0.172, y = 0.1825, label = "Procrustes analysis:", size = 3) +
  annotate(geom = "text", x = 0.172, y = 0.165, label = "italic(t[0]) == 0.9175 *','~italic(p) < 0.0001", parse = TRUE, size = 3) +
  scale_color_manual(values = flPal) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_shape_manual(values = c(24, 25), name = "Depth zone") +
  guides(color = "none", fill = guide_legend(override.aes = list(shape = 15, color = flPal, size = 3), ncol =2, order = 1), shape = guide_legend(ncol = 1)) +
  theme_bw()

sintSymGenProcPlot = sintSymGenProcPlotA +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(color = "black", size = 0.75, fill = NA),
        axis.title = element_blank(),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(color = "black", size = 8),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.box = "horizontal",
        legend.key = element_blank(),
        legend.background = element_blank(),
        legend.title = element_text(color = "black", size = 8),
        legend.text = element_text(color = "black", size = 8)
        )

sintSymGenProcPlot


S. intersepta genetic distance vs Symbiodiniaceae B-C distance

its2DistA = its2Profs %>% arrange(Sample)
its2Dist = its2DistA[c(6:ncol(its2Profs))] %>% decostand("normalize") %>% vegdist(method = "bray")
its2PCoA = cmdscale(its2Dist, eig = TRUE, x.ret = TRUE)

sintIBS = read.table("../data/snps/sintNoClones.ibsMat")[-131,-131] %>% as.matrix() %>% as.dist(diag = FALSE)
sintPCoA = cmdscale(sintIBS, eig = TRUE, x.ret = TRUE)

set.seed(981) 
its2IBSProcrustes = protest(X = sintPCoA, Y = its2PCoA, permutations = 9999)
its2IBSProcrustes
## 
## Call:
## protest(X = sintPCoA, Y = its2PCoA, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.927 
## Correlation in a symmetric Procrustes rotation: 0.2702 
## Significance:  0.0001 
## 
## Permutation: free
## Number of permutations: 9999
plot(its2IBSProcrustes)

plot(its2IBSProcrustes, kind = 2)

admixpops = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone)
admixpops$popdepth = as.factor(paste(admixpops$pop, admixpops$depth, sep = " "))

clumpp4 = read.table("../data/snps/k/ClumppK4.output", header = FALSE)
clumpp4$V1 = admixpops$sample

sintAdmix = admixpops[-131,] %>% left_join(clumpp4[,c(1, 6:9)], by = c("sample" = "V1"))

admixDist = sintAdmix[c(5:ncol(sintAdmix))] %>% vegdist(method = "euclidean")

admixPCoA = cmdscale(admixDist, eig = TRUE, x.ret = TRUE)

set.seed(981) 
its2AdmixProcrustes = protest(X = admixPCoA, Y = its2PCoA, permutations = 9999)

its2AdmixProcrustes
## 
## Call:
## protest(X = admixPCoA, Y = its2PCoA, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.9359 
## Correlation in a symmetric Procrustes rotation: 0.2531 
## Significance:  0.0001 
## 
## Permutation: free
## Number of permutations: 9999
plot(its2AdmixProcrustes)


Cheking dispersion with PERMDISP

Using vegan::betadisper() to look at multivariate homogeneity of dispersion (PERMDISP) between sites and depths. This is using Bray-Curtis dissimilarity.

alpha = with(its2Profs, tapply(specnumber(its2Profs[, c(6:ncol(its2Profs))]), site, mean))
alpha
##  Riley's Hump Tortugas Bank    Lower Keys    Upper Keys 
##      1.177778      1.290909      1.237288      1.200000
gamma = with(its2Profs, specnumber(its2Profs[, c(6:ncol(its2Profs))], site))
gamma
##  Riley's Hump Tortugas Bank    Lower Keys    Upper Keys 
##            13            11            11             9
gamma/alpha
##  Riley's Hump Tortugas Bank    Lower Keys    Upper Keys 
##     11.037736      8.521127      8.890411      7.500000
set.seed(694) 
its2dispS = betadisper(vegdist(decostand(its2Profs[, c(6:ncol(its2Profs))], "normalize")), its2Profs$site)

set.seed(694) 
permutest(its2dispS, permutations = 9999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
## 
## Response: Distances
##            Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups      3 0.0892 0.029731 0.7719   9999  0.516
## Residuals 215 8.2815 0.038518

No significant effect of Site on betadiversity.

alpha = with(its2Profs, tapply(specnumber(its2Profs[, c(6:ncol(its2Profs))]), depthZone, mean))
alpha
##    Shallow Mesophotic 
##   1.310924   1.130000
gamma = with(its2Profs, specnumber(its2Profs[, c(6:ncol(its2Profs))], depthZone))
gamma
##    Shallow Mesophotic 
##         18         11
gamma/alpha
##    Shallow Mesophotic 
##  13.730769   9.734513
set.seed(694)
its2dispD = betadisper(vegdist(decostand(its2Profs[, c(6:ncol(its2Profs))], "normalize")), its2Profs$depthZone)

its2dispD
## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = vegdist(decostand(its2Profs[, c(6:ncol(its2Profs))],
## "normalize")), group = its2Profs$depthZone)
## 
## No. of Positive Eigenvalues: 28
## No. of Negative Eigenvalues: 27
## 
## Average distance to median:
##    Shallow Mesophotic 
##     0.6525     0.5770 
## 
## Eigenvalues for PCoA axes:
## (Showing 8 of 55 eigenvalues)
##  PCoA1  PCoA2  PCoA3  PCoA4  PCoA5  PCoA6  PCoA7  PCoA8 
## 26.118 18.773 13.961  9.281  6.898  6.102  3.901  3.190
set.seed(694)
its2dispDPerm = permutest(its2dispD, permutations = 9999)

its2dispDPerm
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
## 
## Response: Distances
##            Df Sum Sq  Mean Sq      F N.Perm Pr(>F)   
## Groups      1 0.3104 0.310449 10.815   9999 0.0013 **
## Residuals 217 6.2290 0.028705                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Depth does significantly affect beta diversity.

Running PERMANOVA in R

Now let’s see how different communities are from each other with PERMANOVA. We will utilize the vegan::adonis() function. We will use Bray-Curtis similarity for our distance matrix and run a total 0f 9,999 permutations, and test the effects of Site, Depth, and the interaction between Site and Depth.

set.seed(694)
its2Adonis = adonis2(decostand(its2Profs[, c(6:ncol(its2Profs))], "normalize") ~ site * depthZone, data = its2Profs, permutations = 9999, method = "bray")

its2Adonis 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = decostand(its2Profs[, c(6:ncol(its2Profs))], "normalize") ~ site * depthZone, data = its2Profs, permutations = 9999, method = "bray")
##                 Df SumOfSqs      R2       F Pr(>F)    
## site             3    8.455 0.09229  8.4735 0.0001 ***
## depthZone        1    5.368 0.05860 16.1410 0.0001 ***
## site:depthZone   3    7.612 0.08309  7.6290 0.0001 ***
## Residual       211   70.178 0.76602                   
## Total          218   91.613 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(694)
its2PWAdonis = pairwise.adonis(decostand(its2Profs[,c(6:ncol(its2Profs))], "normalize"), factors = its2Profs$site, sim.method = "bray", p.adjust.m = "fdr", perm = 9999)

its2PWAdonis
##                           pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted
## 1 Riley's Hump vs Tortugas Bank  1 3.0054485  7.606522 0.07202701  0.0001    0.00015
## 2    Riley's Hump vs Lower Keys  1 4.5770266 11.695628 0.10286788  0.0001    0.00015
## 3    Riley's Hump vs Upper Keys  1 3.4373577  9.258215 0.08247250  0.0001    0.00015
## 4   Tortugas Bank vs Lower Keys  1 0.9811441  2.446466 0.02137651  0.0334    0.03340
## 5   Tortugas Bank vs Upper Keys  1 1.6940679  4.427003 0.03770004  0.0014    0.00168
## 6      Lower Keys vs Upper Keys  1 3.4239083  9.014882 0.07153823  0.0001    0.00015
##   sig
## 1  **
## 2  **
## 3  **
## 4   .
## 5   *
## 6  **
its2profsRxD = paste(its2Profs$site, its2Profs$depthZone, sep = " ")

set.seed(694)
its2PWAdonis2 = pairwise.adonis(decostand(its2Profs[, c(6:ncol(its2Profs))], "normalize"), factors = its2profsRxD, sim.method = "bray", p.adjust.m = "fdr", perm = 9999)

its2PWAdonis2
##                                                  pairs Df SumsOfSqs   F.Model
## 1      Riley's Hump Shallow vs Riley's Hump Mesophotic  1 0.6958548  1.851042
## 2        Riley's Hump Shallow vs Tortugas Bank Shallow  1 2.8030339  7.429995
## 3     Riley's Hump Shallow vs Tortugas Bank Mesophotic  1 3.0884300  8.640125
## 4           Riley's Hump Shallow vs Lower Keys Shallow  1 5.0368955 15.488866
## 5        Riley's Hump Shallow vs Lower Keys Mesophotic  1 4.6194798 13.467087
## 6           Riley's Hump Shallow vs Upper Keys Shallow  1 4.7577989 13.584574
## 7        Riley's Hump Shallow vs Upper Keys Mesophotic  1 2.2928012  6.845341
## 8     Riley's Hump Mesophotic vs Tortugas Bank Shallow  1 1.7078494  4.482170
## 9  Riley's Hump Mesophotic vs Tortugas Bank Mesophotic  1 1.1310176  3.195846
## 10       Riley's Hump Mesophotic vs Lower Keys Shallow  1 3.2155336 10.357235
## 11    Riley's Hump Mesophotic vs Lower Keys Mesophotic  1 1.8698060  5.584033
## 12       Riley's Hump Mesophotic vs Upper Keys Shallow  1 2.4146387  7.007461
## 13    Riley's Hump Mesophotic vs Upper Keys Mesophotic  1 0.4347921  1.342139
## 14   Tortugas Bank Shallow vs Tortugas Bank Mesophotic  1 2.6960612  7.456037
## 15         Tortugas Bank Shallow vs Lower Keys Shallow  1 1.9880105  6.041735
## 16      Tortugas Bank Shallow vs Lower Keys Mesophotic  1 3.3741943  9.729362
## 17         Tortugas Bank Shallow vs Upper Keys Shallow  1 3.5069150  9.905960
## 18      Tortugas Bank Shallow vs Upper Keys Mesophotic  1 3.2100015  9.476618
## 19      Tortugas Bank Mesophotic vs Lower Keys Shallow  1 5.1052811 16.781413
## 20   Tortugas Bank Mesophotic vs Lower Keys Mesophotic  1 0.3865517  1.192600
## 21      Tortugas Bank Mesophotic vs Upper Keys Shallow  1 1.2421221  3.741093
## 22   Tortugas Bank Mesophotic vs Upper Keys Mesophotic  1 1.2155556  3.855401
## 23         Lower Keys Shallow vs Lower Keys Mesophotic  1 6.2867594 21.368532
## 24            Lower Keys Shallow vs Upper Keys Shallow  1 6.0655965 20.114834
## 25         Lower Keys Shallow vs Upper Keys Mesophotic  1 6.1026797 21.338933
## 26         Lower Keys Mesophotic vs Upper Keys Shallow  1 1.5730255  4.919063
## 27      Lower Keys Mesophotic vs Upper Keys Mesophotic  1 2.7859051  9.149432
## 28         Upper Keys Shallow vs Upper Keys Mesophotic  1 3.3019234 10.593110
##            R2 p.value   p.adjusted sig
## 1  0.04127087  0.1058 0.1139384615    
## 2  0.11355640  0.0001 0.0002000000  **
## 3  0.14017047  0.0001 0.0002000000  **
## 4  0.21367234  0.0001 0.0002000000  **
## 5  0.18843760  0.0001 0.0002000000  **
## 6  0.18976958  0.0001 0.0002000000  **
## 7  0.10556412  0.0003 0.0004941176  **
## 8  0.09439691  0.0008 0.0011200000   *
## 9  0.07757691  0.0243 0.0272160000   .
## 10 0.19781860  0.0001 0.0002000000  **
## 11 0.11493555  0.0025 0.0033333333   *
## 12 0.14012831  0.0003 0.0004941176  **
## 13 0.03026780  0.2359 0.2446370370    
## 14 0.12332989  0.0002 0.0003733333  **
## 15 0.09583706  0.0008 0.0011200000   *
## 16 0.14365058  0.0001 0.0002000000  **
## 17 0.14587762  0.0001 0.0002000000  **
## 18 0.14044299  0.0001 0.0002000000  **
## 19 0.24398180  0.0001 0.0002000000  **
## 20 0.02200670  0.2878 0.2878000000    
## 21 0.06593269  0.0121 0.0147304348   .
## 22 0.06781063  0.0221 0.0257833333   .
## 23 0.27266725  0.0001 0.0002000000  **
## 24 0.26084260  0.0001 0.0002000000  **
## 25 0.27239244  0.0001 0.0002000000  **
## 26 0.07818081  0.0047 0.0059818182   *
## 27 0.13625479  0.0004 0.0006222222  **
## 28 0.15443402  0.0001 0.0002000000  **
its2PWAdonis2 %>% filter(p.adjusted > 0.05)
##                                               pairs Df SumsOfSqs  F.Model         R2
## 1   Riley's Hump Shallow vs Riley's Hump Mesophotic  1 0.6958548 1.851042 0.04127087
## 2  Riley's Hump Mesophotic vs Upper Keys Mesophotic  1 0.4347921 1.342139 0.03026780
## 3 Tortugas Bank Mesophotic vs Lower Keys Mesophotic  1 0.3865517 1.192600 0.02200670
##   p.value p.adjusted sig
## 1  0.1058  0.1139385    
## 2  0.2359  0.2446370    
## 3  0.2878  0.2878000
its2PWAdonis2Tab = its2PWAdonis2 %>% mutate(pairs = pairs, F.Model = round(F.Model, 3), R2 = round(R2,3), p.adjusted =  round(p.adjusted, 4)) %>% select(-p.value, -sig, -SumsOfSqs, -Df) %>%
  flextable() %>%
  flextable::compose(part = "header", j = "pairs", value = as_paragraph("Comparison")) %>%
  flextable::compose(part = "header", j = "F.Model", value = as_paragraph(as_i("Pseudo-F"))) %>%
  flextable::compose(part = "header", j = "R2", value = as_paragraph(as_i("R2"))) %>%
  flextable::compose(part = "header", j = "p.adjusted", value = as_paragraph("p-value")) %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::fontsize(size = 10, part = "all") %>%
  flextable::bold(part = "header") %>%
  flextable::align(align = "left", part = "all") %>%
  flextable::autofit()

table4 = read_docx()
table4 = body_add_flextable(table4, value = its2PWAdonis2Tab)
print(table4, target = "../tables/table4.docx")



Symbiodiniaceae dbRDA

options(sdmpredictors_datadir="../data/snps/bioOracle", timeout = max(300, getOption("timeout")))

flPal = paletteer_d("vapoRwave::jazzCup")[c(2:5)]

popData = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66, 68, 164, 166, 209, 211),] %>% dplyr::select(Sample, site, depthZone, latDD, longDD, depthM)

its2DistA = its2Profs %>% arrange(Sample)
its2Dist = vegdist(decostand(its2DistA[, c(6:ncol(its2DistA))], "normalize"), method = "bray")
datasets = list_datasets(terrestrial = FALSE, marine = TRUE, freshwater = FALSE)

# Extract present-day data sets
present = list_layers(datasets) %>% dplyr::select(dataset_code, layer_code, name, units, description, contains("cellsize"), version) %>% filter(version == 2.2) %>% filter(layer_code %in% c("BO22_damean", "BO22_parmean", "BO22_ph", "BO22_curvelmax_bdmean", "BO22_salinitymean_bdmean", "BO22_salinitymean_ss", "BO22_curvelmean_ss", "BO22_curvelmean_bdmean", "BO22_dissoxmean_bdmean", "BO22_lightbotmax_bdmean", "BO22_lightbotmean_bdmean", "BO22_nitratemean_bdmean", "BO22_tempmax_bdmean", "BO22_tempmean_bdmean", "BO22_ppmean_ss", "BO22_ppmean_bdmean", "BO22_chlomean_ss", "BO22_chlomean_bdmean"))

envVar = load_layers(present$layer_code)

itsEnvData2 = data.frame(popData, raster::extract(envVar, popData[,5:4]))[-131,]

corData = rcorr(as.matrix(itsEnvData2[,c(6:ncol(itsEnvData2))]), type = "pearson")
corDataFlat = melt(corData$r, value.name = "r")
pDataFlat = melt(corData$P, value.name = "p")
corDataBind = corDataFlat %>% left_join(pDataFlat, by = c("Var1","Var2"))

ggplot(corDataBind) +
  geom_tile(aes(x = Var1, y = Var2, fill = r)) +
  scale_fill_gradient2(low = "#3B9AB2FF", high = "#F21A00FF") +
  geom_text(data = filter(corDataBind, r >= 0.7, p < 0.05),aes(x = Var1, y = Var2, label = round(r, 2))) +
    geom_text(data = filter(corDataBind, r <= -0.7, p < 0.05),aes(x = Var1, y = Var2, label = round(r, 2))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

itsEnvData2 = itsEnvData2 %>% dplyr::select("BO22_curvelmean_ss", "depthM", "BO22_lightbotmean_bdmean", "BO22_tempmean_bdmean", "BO22_ppmean_bdmean", "BO22_dissoxmean_bdmean")
corData2 = cor(itsEnvData2)
corMelt2 = melt(corData2)

ggplot(corMelt2) +
  geom_tile(aes(x = Var1, y = Var2, fill = value)) +
  scale_fill_gradient2(low = "#3B9AB2FF", high = "#F21A00FF") +
  geom_text(data = corMelt2[corMelt2$value >= 0.7,],aes(x = Var1, y = Var2, label = round(value, 2))) +
    geom_text(data = corMelt2[corMelt2$value <= -0.7,],aes(x = Var1, y = Var2, label = round(value, 2))) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90))

vif = diag(solve(cor(itsEnvData2)))
vif
##       BO22_curvelmean_ss                   depthM BO22_lightbotmean_bdmean 
##                 4.023795                 1.693764                 8.317402 
##     BO22_tempmean_bdmean       BO22_ppmean_bdmean   BO22_dissoxmean_bdmean 
##                 6.916539                 8.158537                 6.199858
itsDbrda0 = dbrda(its2Dist ~ 1, data = itsEnvData2)
itsDbrda1 = dbrda(its2Dist ~ ., data = itsEnvData2)
anova(itsDbrda1)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ BO22_curvelmean_ss + depthM + BO22_lightbotmean_bdmean + BO22_tempmean_bdmean + BO22_ppmean_bdmean + BO22_dissoxmean_bdmean, data = itsEnvData2)
##           Df SumOfSqs      F Pr(>F)    
## Model      6   14.774 6.7937  0.001 ***
## Residual 212   76.839                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(itsDbrda1)
## $r.squared
## [1] 0.1612662
## 
## $adj.r.squared
## [1] 0.1375284
set.seed(092)
bestItsDbrda <- ordiR2step(itsDbrda0, itsDbrda1)
## Step: R2.adj= 0 
## Call: its2Dist ~ 1 
##  
##                            R2.adjusted
## <All variables>             0.13752844
## + depthM                    0.07556852
## + BO22_curvelmean_ss        0.03225721
## + BO22_dissoxmean_bdmean    0.02317990
## + BO22_tempmean_bdmean      0.02022712
## + BO22_ppmean_bdmean        0.01667414
## + BO22_lightbotmean_bdmean  0.01473131
## <none>                      0.00000000
## 
##          Df    AIC      F Pr(>F)   
## + depthM  1 974.13 18.821  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.07556852 
## Call: its2Dist ~ depthM 
##  
##                            R2.adjusted
## <All variables>             0.13752844
## + BO22_curvelmean_ss        0.10492630
## + BO22_lightbotmean_bdmean  0.09605059
## + BO22_ppmean_bdmean        0.09267258
## + BO22_tempmean_bdmean      0.09002396
## + BO22_dissoxmean_bdmean    0.08748718
## <none>                      0.07556852
## 
##                      Df    AIC      F Pr(>F)   
## + BO22_curvelmean_ss  1 968.05 8.1174  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1049263 
## Call: its2Dist ~ depthM + BO22_curvelmean_ss 
##  
##                            R2.adjusted
## <All variables>              0.1375284
## + BO22_lightbotmean_bdmean   0.1255916
## + BO22_tempmean_bdmean       0.1210367
## + BO22_ppmean_bdmean         0.1189578
## + BO22_dissoxmean_bdmean     0.1175131
## <none>                       0.1049263
## 
##                            Df    AIC      F Pr(>F)   
## + BO22_lightbotmean_bdmean  1 963.92 6.1048  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1255916 
## Call: its2Dist ~ depthM + BO22_curvelmean_ss + BO22_lightbotmean_bdmean 
##  
##                          R2.adjusted
## <All variables>            0.1375284
## + BO22_dissoxmean_bdmean   0.1329506
## + BO22_tempmean_bdmean     0.1310355
## + BO22_ppmean_bdmean       0.1272419
## <none>                     0.1255916
## 
##                          Df    AIC      F Pr(>F)   
## + BO22_dissoxmean_bdmean  1 963.05 2.8248   0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1329506 
## Call: its2Dist ~ depthM + BO22_curvelmean_ss + BO22_lightbotmean_bdmean +      BO22_dissoxmean_bdmean 
##  
##                        R2.adjusted
## <All variables>          0.1375284
## + BO22_tempmean_bdmean   0.1347526
## <none>                   0.1329506
## + BO22_ppmean_bdmean     0.1323477
## 
##                        Df    AIC      F Pr(>F)
## + BO22_tempmean_bdmean  1 963.57 1.4457  0.158
bestItsDbrda$anova
##                              R2.adj Df    AIC       F Pr(>F)   
## + depthM                   0.075569  1 974.13 18.8206  0.002 **
## + BO22_curvelmean_ss       0.104926  1 968.05  8.1174  0.002 **
## + BO22_lightbotmean_bdmean 0.125592  1 963.92  6.1048  0.002 **
## + BO22_dissoxmean_bdmean   0.132951  1 963.05  2.8248  0.010 **
## <All variables>            0.137528                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now create best model and look at variance partitioning and significance with ANOVA

its2Model = itsEnvData2 %>% dplyr::select("Depth" = depthM, "Current" = BO22_curvelmean_ss, "DO" = BO22_dissoxmean_bdmean, "Light" = BO22_lightbotmean_bdmean)

its2Dbrda = dbrda(its2Dist ~ Depth + Current + DO + Light, its2Model)

itsDbrdaVarPart = varpart(its2Dist, ~ depthM, ~ BO22_curvelmean_ss, ~ BO22_dissoxmean_bdmean, ~ BO22_lightbotmean_bdmean, data = itsEnvData2)

itsDbrdaDepth = dbrda(its2Dist ~ Depth, its2Model)
itsDbrdaCur = dbrda(its2Dist ~ Current, its2Model)
itsDbrdaDO = dbrda(its2Dist ~ DO, its2Model)
itsDbrdaLight = dbrda(its2Dist ~ Light, its2Model)

plot(itsDbrdaVarPart, Xnames = c("Depth", "Current", "DO", "Light"), bg = c(4:7), digits = 2)

itsDbrdaVarPart
## 
## Partition of squared Bray distance in dbRDA 
## 
## Call: varpart(Y = its2Dist, X = ~depthM, ~BO22_curvelmean_ss,
## ~BO22_dissoxmean_bdmean, ~BO22_lightbotmean_bdmean, data = itsEnvData2)
## 
## Explanatory tables:
## X1:  ~depthM
## X2:  ~BO22_curvelmean_ss
## X3:  ~BO22_dissoxmean_bdmean
## X4:  ~BO22_lightbotmean_bdmean 
## 
## No. of explanatory tables: 4 
## Total variation (SS): 91.613 
## No. of observations: 219 
## 
## Partition table:
##                             Df R.square Adj.R.square Testable
## [aeghklno] = X1              1  0.07981      0.07557     TRUE
## [befiklmo] = X2              1  0.03670      0.03226     TRUE
## [cfgjlmno] = X3              1  0.02766      0.02318     TRUE
## [dhijkmno] = X4              1  0.01925      0.01473     TRUE
## [abefghiklmno] = X1+X2       2  0.11314      0.10493     TRUE
## [acefghjklmno] = X1+X3       2  0.09586      0.08749     TRUE
## [adeghijklmno] = X1+X4       2  0.10434      0.09605     TRUE
## [bcefgijklmno] = X2+X3       2  0.06500      0.05635     TRUE
## [bdefhijklmno] = X2+X4       2  0.05629      0.04755     TRUE
## [cdfghijklmno] = X3+X4       2  0.04408      0.03523     TRUE
## [abcefghijklmno] = X1+X2+X3  3  0.12966      0.11751     TRUE
## [abdefghijklmno] = X1+X2+X4  3  0.13762      0.12559     TRUE
## [acdefghijklmno] = X1+X3+X4  3  0.11504      0.10269     TRUE
## [bcdefghijklmno] = X2+X3+X4  3  0.08217      0.06937     TRUE
## [abcdefghijklmno] = All      4  0.14886      0.13295     TRUE
## Individual fractions                                         
## [a] = X1 | X2+X3+X4          1               0.06358     TRUE
## [b] = X2 | X1+X3+X4          1               0.03026     TRUE
## [c] = X3 | X1+X2+X4          1               0.00736     TRUE
## [d] = X4 | X1+X2+X3          1               0.01544     TRUE
## [e]                          0               0.00388    FALSE
## [f]                          0              -0.00072    FALSE
## [g]                          0               0.01445    FALSE
## [h]                          0              -0.00242    FALSE
## [i]                          0              -0.00023    FALSE
## [j]                          0               0.00523    FALSE
## [k]                          0              -0.00074    FALSE
## [l]                          0              -0.00060    FALSE
## [m]                          0               0.00005    FALSE
## [n]                          0              -0.00295    FALSE
## [o]                          0               0.00036    FALSE
## [p] = Residuals              0               0.86705    FALSE
## Controlling 2 tables X                                       
## [ae] = X1 | X3+X4            1               0.06746     TRUE
## [ag] = X1 | X2+X4            1               0.07804     TRUE
## [ah] = X1 | X2+X3            1               0.06117     TRUE
## [be] = X2 | X3+X4            1               0.03414     TRUE
## [bf] = X2 | X1+X4            1               0.02954     TRUE
## [bi] = X2 | X1+X3            1               0.03003     TRUE
## [cf] = X3 | X1+X4            1               0.00664     TRUE
## [cg] = X3 | X2+X4            1               0.02181     TRUE
## [cj] = X3 | X1+X2            1               0.01259     TRUE
## [dh] = X4 | X2+X3            1               0.01302     TRUE
## [di] = X4 | X1+X3            1               0.01520     TRUE
## [dj] = X4 | X1+X2            1               0.02067     TRUE
## Controlling 1 table X                                        
## [aghn] = X1 | X2             1               0.07267     TRUE
## [aehk] = X1 | X3             1               0.06431     TRUE
## [aegl] = X1 | X4             1               0.08132     TRUE
## [bfim] = X2 | X1             1               0.02936     TRUE
## [beik] = X2 | X3             1               0.03317     TRUE
## [befl] = X2 | X4             1               0.03282     TRUE
## [cfjm] = X3 | X1             1               0.01192     TRUE
## [cgjn] = X3 | X2             1               0.02409     TRUE
## [cfgl] = X3 | X4             1               0.02050     TRUE
## [dijm] = X4 | X1             1               0.02048     TRUE
## [dhjn] = X4 | X2             1               0.01530     TRUE
## [dhik] = X4 | X3             1               0.01205     TRUE
## ---
## Use function 'dbrda' to test significance of fractions of interest
set.seed(004)
anova(its2Dbrda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ Depth + Current + DO + Light, data = its2Model)
##           Df SumOfSqs      F Pr(>F)    
## Model      4   13.638 9.3569  0.001 ***
## Residual 214   77.976                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(003)
anova(itsDbrdaDepth)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ Depth, data = its2Model)
##           Df SumOfSqs      F Pr(>F)    
## Model      1    7.312 18.821  0.001 ***
## Residual 217   84.302                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(002)
anova(itsDbrdaCur)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ Current, data = its2Model)
##           Df SumOfSqs      F Pr(>F)    
## Model      1    3.362 8.2665  0.001 ***
## Residual 217   88.251                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(001)
anova(itsDbrdaDO)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ DO, data = its2Model)
##           Df SumOfSqs      F Pr(>F)    
## Model      1    2.534 6.1731  0.001 ***
## Residual 217   89.079                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(000)
anova(itsDbrdaLight)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = its2Dist ~ Light, data = its2Model)
##           Df SumOfSqs      F Pr(>F)    
## Model      1    1.764 4.2594  0.001 ***
## Residual 217   89.850                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prepare data and plot dbRDA

itsRdaVar = round(its2Dbrda$CA$eig/sum(its2Dbrda$CA$eig)*100, 1)
head(itsRdaVar)
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 
## 31.3 19.2 14.6  9.6  8.1  5.9
itsRdaVarFit = round(its2Dbrda$CCA$eig/sum(its2Dbrda$CCA$eig)*100, 1)
head(itsRdaVarFit)
## dbRDA1 dbRDA2 dbRDA3 dbRDA4 
##   59.2   21.1   15.0    4.7
sintI2P = read.csv("../data/stephanocoeniaMetaData.csv")[-c(66,68,164,166,133,209,211),] %>% dplyr::select("sample" = tubeID, "pop" = site, "depth" = depthZone)

sintI2P$popdepth = paste(sintI2P$pop, sintI2P$depth)

its2RdaPoints = as.data.frame(scores(its2Dbrda, choices=c(1:4))$sites)
its2RdaPoints$sample = sintI2P$sample
head(its2RdaPoints)
##        dbRDA1     dbRDA2      dbRDA3     dbRDA4 sample
## 1 -1.31749961  2.7339073  0.19012497  2.0456033 SFK001
## 2 -0.08748044 -0.9340035  1.27584673 -2.6282432 SFK002
## 3 -0.09302170 -0.9302237  1.30967725 -2.5170507 SFK003
## 4 -1.20478994  2.4966246  0.12686050  1.4064890 SFK004
## 5 -0.82873686 -2.7158800  0.09158796  0.5985444 SFK005
## 6 -0.37542179  1.5084117 -0.72284403 -1.6586371 SFK006
its2DbrdaData1 = sintI2P %>% left_join(its2RdaPoints)
## Joining, by = "sample"
head(its2DbrdaData1)
##   sample           pop      depth                 popdepth      dbRDA1     dbRDA2
## 1 SFK001 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -1.31749961  2.7339073
## 2 SFK002 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -0.08748044 -0.9340035
## 3 SFK003 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -0.09302170 -0.9302237
## 4 SFK004 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -1.20478994  2.4966246
## 5 SFK005 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -0.82873686 -2.7158800
## 6 SFK006 Tortugas Bank Mesophotic Tortugas Bank Mesophotic -0.37542179  1.5084117
##        dbRDA3     dbRDA4
## 1  0.19012497  2.0456033
## 2  1.27584673 -2.6282432
## 3  1.30967725 -2.5170507
## 4  0.12686050  1.4064890
## 5  0.09158796  0.5985444
## 6 -0.72284403 -1.6586371
tail(its2DbrdaData1)
##     sample        pop      depth              popdepth     dbRDA1     dbRDA2
## 214 SFK215 Upper Keys Mesophotic Upper Keys Mesophotic -1.3174996  2.7339073
## 215 SFK216 Upper Keys Mesophotic Upper Keys Mesophotic -1.3174996  2.7339073
## 216 SFK217 Upper Keys    Shallow    Upper Keys Shallow -0.8287369 -2.7158800
## 217 SFK218 Upper Keys    Shallow    Upper Keys Shallow  0.5921418 -0.1975787
## 218 SFK219 Upper Keys    Shallow    Upper Keys Shallow -0.8287369 -2.7158800
## 219 SFK220 Upper Keys    Shallow    Upper Keys Shallow -0.0930217 -0.9302237
##          dbRDA3     dbRDA4
## 214  0.19012497  2.0456033
## 215  0.19012497  2.0456033
## 216  0.09158796  0.5985444
## 217 -5.38330662  1.7026079
## 218  0.09158796  0.5985444
## 219  1.30967725 -2.5170507
its2EnvLoad = as.data.frame(its2Dbrda$CCA$biplot)
its2EnvLoad$var = row.names(its2EnvLoad) 
its2EnvLoad$var = its2EnvLoad$var %>% gsub(pattern = "`", replacement = "", x = its2EnvLoad$var)

its2DbrdaData = merge(its2DbrdaData1, aggregate(cbind(mean.x = dbRDA1, mean.y2 = dbRDA2, mean.y3 = dbRDA3, mean.y4 = dbRDA4)~popdepth, its2DbrdaData1, mean), by="popdepth")

its2DbrdaData$depth = factor(its2DbrdaData$depth)
its2DbrdaData$depth = factor(its2DbrdaData$depth, levels(its2DbrdaData$depth)[c(2,1)])
its2DbrdaData$pop = factor(its2DbrdaData$pop)
its2DbrdaData$pop = factor(its2DbrdaData$pop, levels(its2DbrdaData$pop)[c(4, 1, 3, 2)])

head(its2DbrdaData)
##                popdepth sample        pop      depth      dbRDA1     dbRDA2     dbRDA3
## 1 Lower Keys Mesophotic SFK101 Lower Keys Mesophotic -0.09945113 -0.9528654 1.21166214
## 2 Lower Keys Mesophotic SFK102 Lower Keys Mesophotic -0.82873686 -2.7158800 0.09158796
## 3 Lower Keys Mesophotic SFK103 Lower Keys Mesophotic -1.31749961  2.7339073 0.19012497
## 4 Lower Keys Mesophotic SFK104 Lower Keys Mesophotic -0.09302170 -0.9302237 1.30967725
## 5 Lower Keys Mesophotic SFK105 Lower Keys Mesophotic -0.83535932 -2.7450672 0.04917306
## 6 Lower Keys Mesophotic SFK106 Lower Keys Mesophotic -0.82873686 -2.7158800 0.09158796
##       dbRDA4     mean.x mean.y2   mean.y3   mean.y4
## 1 -2.9335212 -0.4011398 -1.5888 0.4928093 -1.347814
## 2  0.5985444 -0.4011398 -1.5888 0.4928093 -1.347814
## 3  2.0456033 -0.4011398 -1.5888 0.4928093 -1.347814
## 4 -2.5170507 -0.4011398 -1.5888 0.4928093 -1.347814
## 5  0.2685841 -0.4011398 -1.5888 0.4928093 -1.347814
## 6  0.5985444 -0.4011398 -1.5888 0.4928093 -1.347814
its2DbrdaPlotA = ggplot() +
  geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
  geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
  geom_segment(data = its2EnvLoad, aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2), color = "gray35", arrow = arrow(length = unit(0.15, "cm"), type = "open"), size = 0.65) +
  geom_point(data = its2DbrdaData, aes(x = dbRDA1, y = dbRDA2, fill = pop, shape = depth, color = pop), size = 2, alpha = 0.5, show.legend = FALSE, position = position_jitter(seed = 1, width = 0.075, height = 0.075)) +
  scale_shape_manual(values = c(24, 25), name = "Depth Zone") +
  geom_point(data = its2DbrdaData, aes(x = mean.x, y = mean.y2, fill = pop, shape = depth), size = 3, color = "black") + #population centroids indicated by triangles
  geom_text(data = its2EnvLoad[1,], aes(x = dbRDA1-0.2, y = dbRDA2-0.02, label = var), color = "gray35", size = 3) +
  geom_text(data = its2EnvLoad[2,], aes(x = dbRDA1+0.1, y = dbRDA2+0.2, label = var), color = "gray35", size = 3) +
  geom_text(data = its2EnvLoad[3,], aes(x = dbRDA1+0.15, y = dbRDA2+0.1, label = var), color = "gray35", size = 3) +
  geom_text(data = its2EnvLoad[4,], aes(x = dbRDA1, y = dbRDA2-0.1, label = var), color = "gray35", size = 3) +
  scale_fill_manual(values = flPal, name = "Site") +
  scale_color_manual(values = flPal, name = "Site") +
  labs(x = paste0("dbRDA 1 (", itsRdaVarFit[1], "% [",itsRdaVar[1], "%])"), y = paste0("dbRDA 2 (", itsRdaVarFit[2], "% [",itsRdaVar[2], "%])")) +
  guides(shape = guide_legend(override.aes = list(size = 3, stroke = 0.5, alpha = NA), order = 2), fill = guide_legend(override.aes = list(shape = 22, size = 4, color = flPal, alpha = NA), order = 1))+
  theme_bw()

its2DbrdaPlot = its2DbrdaPlotA +
  theme(axis.title.x = element_text(color = "black", size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank(),
        axis.title.y = element_text(color = "black", size = 10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        legend.position = "right",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8),
        panel.border = element_rect(color = "black", size = 1),
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

its2DbrdaPlot

ggsave("../figures/figure8.png", plot = its2DbrdaPlot, height = 8, width = 14, units = "cm", dpi = 300)
ggsave("../figures/figure8.svg", plot = its2DbrdaPlot, height = 8, width = 14, units = "cm", dpi = 300)